diff --git a/applications/test/extendedStencil/Test-ExtendedStencil2.C b/applications/test/extendedStencil/Test-ExtendedStencil2.C new file mode 100644 index 0000000000000000000000000000000000000000..eda8f986793c781bfc86f61a26667b7b03642fd8 --- /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 <http://www.gnu.org/licenses/>. + +Application + testExtendedStencil2 + +Description + Test app for determining extended cell-to-cell stencil. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "fvMesh.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "Time.H" +#include "OFstream.H" +#include "meshTools.H" + +#include "CFCCellToCellStencil.H" + + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void writeStencilOBJ +( + const fileName& fName, + const point& fc, + const List<point>& stencilCc +) +{ + OFstream str(fName); + label vertI = 0; + + meshTools::writeOBJ(str, fc); + vertI++; + + forAll(stencilCc, i) + { + meshTools::writeOBJ(str, stencilCc[i]); + vertI++; + str << "l 1 " << vertI << nl; + } +} + + +// Stats +void writeStencilStats(const labelListList& stencil) +{ + label sumSize = 0; + label nSum = 0; + label minSize = labelMax; + label maxSize = labelMin; + + forAll(stencil, i) + { + const labelList& sCells = stencil[i]; + + if (sCells.size() > 0) + { + sumSize += sCells.size(); + nSum++; + minSize = min(minSize, sCells.size()); + maxSize = max(maxSize, sCells.size()); + } + } + reduce(sumSize, sumOp<label>()); + reduce(nSum, sumOp<label>()); + sumSize /= nSum; + + reduce(minSize, minOp<label>()); + reduce(maxSize, maxOp<label>()); + + Info<< "Stencil size :" << nl + << " average : " << sumSize << nl + << " min : " << minSize << nl + << " max : " << maxSize << nl + << endl; +} + + +// Main program: + +int main(int argc, char *argv[]) +{ +# include "addTimeOptions.H" +# include "setRootCase.H" +# include "createTime.H" + + // Get times list + instantList Times = runTime.times(); +# include "checkTimeOptions.H" + runTime.setTime(Times[startTime], startTime); +# include "createMesh.H" + + + //---- CELL CENTRED STENCIL ----- + + // Centred, cell stencil + // ~~~~~~~~~~~~~~~~~~~~~ + + { + // Full stencil. This is per local cell a set of global indices, either + // into cells or also boundary faces. + CFCCellToCellStencil stencil(mesh); + + // Construct exchange map and renumber stencil + List<Map<label> > compactMap(Pstream::nProcs()); + mapDistribute map + ( + stencil.globalNumbering(), + stencil, + compactMap + ); + + + // Print some stats + Info<< "cellFaceCell:" << endl; + writeStencilStats(stencil); + + + // Collect the data in stencil form + List<List<vector> > stencilPoints; + { + const volVectorField& fld = mesh.C(); + + // 1. Construct cell data in compact addressing + List<point> compactFld(map.constructSize(), pTraits<point>::zero); + + // Insert my internal values + forAll(fld, cellI) + { + compactFld[cellI] = fld[cellI]; + } + // Insert my boundary values + label nCompact = fld.size(); + forAll(fld.boundaryField(), patchI) + { + const fvPatchField<vector>& pfld = fld.boundaryField()[patchI]; + + forAll(pfld, i) + { + compactFld[nCompact++] = pfld[i]; + } + } + + // Do all swapping + map.distribute(compactFld); + + // 2. Pull to stencil + stencilPoints.setSize(stencil.size()); + + forAll(stencil, cellI) + { + const labelList& compactCells = stencil[cellI]; + + stencilPoints[cellI].setSize(compactCells.size()); + + forAll(compactCells, i) + { + stencilPoints[cellI][i] = compactFld[compactCells[i]]; + } + } + } + + + forAll(stencilPoints, cellI) + { + writeStencilOBJ + ( + runTime.path()/"centredCell" + Foam::name(cellI) + ".obj", + mesh.cellCentres()[cellI], + stencilPoints[cellI] + ); + } + } + + Pout<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C b/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C index 14959901f7921e634ceabab35fa23f2321c9c691..98ead1f4bbe808d2676ce1e352f24e41e69e3887 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<word>& selectedFields_; instant ti_; instant ti1_; + const interpolationWeights& interpolator_; + const wordList& timeNames_; int divisions_; public: @@ -60,6 +64,8 @@ public: const HashSet<word>& selectedFields, const instant& ti, const instant& ti1, + const interpolationWeights& interpolator, + const wordList& timeNames, int divisions ) : @@ -69,6 +75,8 @@ public: selectedFields_(selectedFields), ti_(ti), ti1_(ti1), + interpolator_(interpolator), + timeNames_(timeNames), divisions_(divisions) {} @@ -98,34 +106,6 @@ void fieldInterpolator::interpolate() { Info<< " " << fieldIter()->name() << '('; - GeoFieldType fieldi - ( - IOobject - ( - fieldIter()->name(), - ti_.name(), - fieldIter()->db(), - IOobject::MUST_READ, - IOobject::NO_WRITE, - false - ), - mesh_ - ); - - GeoFieldType fieldi1 - ( - IOobject - ( - fieldIter()->name(), - ti1_.name(), - fieldIter()->db(), - IOobject::MUST_READ, - IOobject::NO_WRITE, - false - ), - mesh_ - ); - scalar deltaT = (ti1_.value() - ti_.value())/(divisions_ + 1); for (int j=0; j<divisions_; j++) @@ -141,20 +121,51 @@ void fieldInterpolator::interpolate() Info<< " "; } - scalar lambda = scalar(j + 1)/scalar(divisions_ + 1); + // Calculate times to read and weights + labelList indices; + scalarField weights; + interpolator_.valueWeights + ( + runTime_.value(), + indices, + weights + ); + + const wordList selectedTimeNames + ( + UIndirectList<word>(timeNames_, indices)() + ); + + //Info<< "For time " << runTime_.value() + // << " need times " << selectedTimeNames + // << " need weights " << weights << endl; + + + // Read on the objectRegistry all the required fields + ReadFields<GeoFieldType> + ( + fieldIter()->name(), + mesh_, + selectedTimeNames + ); GeoFieldType fieldj ( - IOobject + uniformInterpolate<GeoFieldType> ( + IOobject + ( + fieldIter()->name(), + runTime_.timeName(), + fieldIter()->db(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), fieldIter()->name(), - timej.name(), - fieldIter()->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - (1.0 - lambda)*fieldi + lambda*fieldi1 + selectedTimeNames, + weights + ) ); fieldj.write(); @@ -188,6 +199,12 @@ int main(int argc, char *argv[]) "integer", "specify number of temporal sub-divisions to create (default = 1)." ); + argList::addOption + ( + "interpolationType", + "word", + "specify type of interpolation (linear or spline)" + ); #include "setRootCase.H" #include "createTime.H" @@ -198,15 +215,51 @@ int main(int argc, char *argv[]) { args.optionLookup("fields")() >> selectedFields; } + if (selectedFields.size()) + { + Info<< "Interpolating fields " << selectedFields << nl << endl; + } + else + { + Info<< "Interpolating all fields" << nl << endl; + } + int divisions = 1; if (args.optionFound("divisions")) { args.optionLookup("divisions")() >> divisions; } + Info<< "Using " << divisions << " per time interval" << nl << endl; + + + const word interpolationType = args.optionLookupOrDefault<word> + ( + "interpolationType", + "linear" + ); + Info<< "Using interpolation " << interpolationType << nl << endl; + instantList timeDirs = timeSelector::select0(runTime, args); + scalarField timeVals(timeDirs.size()); + wordList timeNames(timeDirs.size()); + forAll(timeDirs, i) + { + timeVals[i] = timeDirs[i].value(); + timeNames[i] = timeDirs[i].name(); + } + autoPtr<interpolationWeights> interpolatorPtr + ( + interpolationWeights::New + ( + interpolationType, + timeVals + ) + ); + + #include "createMesh.H" Info<< "Interpolating fields for times:" << endl; @@ -226,6 +279,8 @@ int main(int argc, char *argv[]) selectedFields, timeDirs[timei], timeDirs[timei+1], + interpolatorPtr(), + timeNames, divisions ); diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index c58cf1afc292006cb18bc77bdd9e82540d2663c7..948ec4dbe7a3e08cfeea1fb955c645e0a986a582 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/containers/Lists/UIndirectList/UIndirectList.H b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H index 33e929baccdfe64220b30295a9ce7c7cd28621c3..e33acfa4643c3a01029c5300ba41759f6e74790e 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 diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.C b/src/OpenFOAM/db/objectRegistry/objectRegistry.C index 6fcaed2e2eab3df27a80a878edccd0d3e469ed0b..09307c3a5ce20194f5822e604829133d65748826 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<objectRegistry>(name)) + { + objectRegistry* fieldsCachePtr = new objectRegistry + ( + IOobject + ( + name, + time().constant(), + *this, + IOobject::NO_READ, + IOobject::NO_WRITE + ) + ); + fieldsCachePtr->store(); + } return lookupObject<objectRegistry>(name); } diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.H b/src/OpenFOAM/db/objectRegistry/objectRegistry.H index 34ddc8fc348946e1a2869d527ad27cd80dbd2582..48cabcc1cfa79c63242306fb2654e3099cecaf59 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<class Type> wordList names() const; - //- Lookup and return a const sub-objectRegistry - const objectRegistry& subRegistry(const word& name) const; + //- Lookup and return a const sub-objectRegistry. Optionally create + // it if it does not exist. + const objectRegistry& subRegistry + ( + const word& name, + const bool forceCreate = false + ) const; //- Lookup and return all objects of the given Type template<class Type> diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C new file mode 100644 index 0000000000000000000000000000000000000000..968f5f4c0f0529d5dce4d573359b0b524c53761c --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class GeoField> +Foam::tmp<GeoField> Foam::uniformInterpolate +( + const HashPtrTable<GeoField, label, Hash<label> >& fields, + const labelList& indices, + const scalarField& weights +) +{ + const GeoField& field0 = *(*fields.begin()); + + // Interpolate + tmp<GeoField> tfld + ( + new GeoField + ( + IOobject + ( + "uniformInterpolate(" + field0.name() + ')', + field0.time().timeName(), + field0.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + weights[0]*(*fields[indices[0]]) + ) + ); + GeoField& fld = tfld(); + + for (label i = 1; i < indices.size(); ++i) + { + fld += weights[i]*(*fields[indices[i]]); + } + + return tfld; +} + + +template<class GeoField> +Foam::tmp<GeoField> Foam::uniformInterpolate +( + const IOobject& fieldIO, + const word& fieldName, + const wordList& times, + const scalarField& weights, + const objectRegistry& fieldsCache +) +{ + // Look up the first field + const objectRegistry& time0Fields = fieldsCache.lookupObject + < + const objectRegistry + > + ( + times[0] + ); + const GeoField& field0 = time0Fields.lookupObject + < + const GeoField + > + ( + fieldName + ); + + + // Interpolate + tmp<GeoField> tfld(new GeoField(fieldIO, weights[0]*field0)); + GeoField& fld = tfld(); + + for (label i = 1; i < times.size(); ++i) + { + const objectRegistry& timeIFields = fieldsCache.lookupObject + < + const objectRegistry + > + ( + times[i] + ); + const GeoField& fieldI = timeIFields.lookupObject + < + const GeoField + > + ( + fieldName + ); + + fld += weights[i]*fieldI; + } + + return tfld; +} + + +template<class GeoField> +Foam::tmp<GeoField> Foam::uniformInterpolate +( + const IOobject& fieldIO, + const word& fieldName, + const wordList& times, + const scalarField& weights, + const word& registryName +) +{ + return uniformInterpolate<GeoField> + ( + fieldIO, + fieldName, + times, + weights, + fieldIO.db().subRegistry(registryName, true) + ); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H new file mode 100644 index 0000000000000000000000000000000000000000..f196f71c192e64e4bae65f268a2e57396b048d42 --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "GeometricField.H" +#include "HashPtrTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * Global functions * * * * * * * * * * * * * // + +//- Interpolate selected fields (given by indices and corresponding weights) +// (boundary type becomes calculated). Fields stored per index. Field gets name +// "uniformInterpolate(" + fld.name() + ')'. +template<class GeoField> +tmp<GeoField> uniformInterpolate +( + const HashPtrTable<GeoField, label, Hash<label> >& fields, + const labelList& indices, + const scalarField& weights +); + +//- Interpolate fields. fieldsCache contains per timeName all loaded fields. +// Resulting field gets properties according to fieldIO +template<class GeoField> +tmp<GeoField> uniformInterpolate +( + const IOobject& fieldIO, + const word& fieldName, + const wordList& times, + const scalarField& weights, + const objectRegistry& fieldsCache +); + +//- Interpolate fields. fieldsCache contains per timeName all loaded fields. +// Resulting field gets properties according to fieldIO +template<class GeoField> +tmp<GeoField> uniformInterpolate +( + const IOobject& fieldIO, + const word& fieldName, + const wordList& times, + const scalarField& weights, + const word& registryName = "fieldsCache" +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "uniformInterpolate.C" +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/fields/ReadFields/ReadFields.C b/src/OpenFOAM/fields/ReadFields/ReadFields.C index 7a17b362dc3a5eaa568081840ffddde5884f452f..a679082b417626c95f1928af2cc1c84ee8dd3b78 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<class GeoField> +void Foam::ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + objectRegistry& fieldsCache +) +{ + // Collect all times that are no longer used + { + HashSet<word> usedTimes(timeNames); + + DynamicList<word> unusedTimes(fieldsCache.size()); + + forAllIter(objectRegistry, fieldsCache, timeIter) + { + const word& tm = timeIter.key(); + if (!usedTimes.found(tm)) + { + unusedTimes.append(tm); + } + } + + //Info<< "Unloading times " << unusedTimes << endl; + + forAll(unusedTimes, i) + { + objectRegistry& timeCache = const_cast<objectRegistry&> + ( + fieldsCache.lookupObject<objectRegistry>(unusedTimes[i]) + ); + fieldsCache.checkOut(timeCache); + } + } + + + // Load any new fields + forAll(timeNames, i) + { + const word& tm = timeNames[i]; + + // Create if not found + if (!fieldsCache.found(tm)) + { + //Info<< "Creating registry for time " << tm << endl; + + // Create objectRegistry if not found + objectRegistry* timeCachePtr = new objectRegistry + ( + IOobject + ( + tm, + tm, + fieldsCache, + IOobject::NO_READ, + IOobject::NO_WRITE + ) + ); + timeCachePtr->store(); + } + + // Obtain cache for current time + const objectRegistry& timeCache = + fieldsCache.lookupObject<objectRegistry> + ( + tm + ); + + // Store field if not found + if (!timeCache.found(fieldName)) + { + //Info<< "Loading field " << fieldName + // << " for time " << tm << endl; + + GeoField loadedFld + ( + IOobject + ( + fieldName, + tm, + mesh.thisDb(), + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ), + mesh + ); + + // Transfer to timeCache (new objectRegistry and store flag) + GeoField* fldPtr = new GeoField + ( + IOobject + ( + fieldName, + tm, + timeCache, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + loadedFld + ); + fldPtr->store(); + } + } +} + + +template<class GeoField> +void Foam::ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + const word& registryName +) +{ + ReadFields<GeoField> + ( + fieldName, + mesh, + timeNames, + const_cast<objectRegistry&> + ( + mesh.thisDb().subRegistry(registryName, true) + ) + ); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/fields/ReadFields/ReadFields.H b/src/OpenFOAM/fields/ReadFields/ReadFields.H index 11ff547935c953a0d304addb781c0f2f204bd0d1..ac4bdb37c740e385ee3ad4b1a14ceadb0459ceca 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<class GeoField> +static void ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + objectRegistry& fieldsCache +); + +//- Helper routine to read GeometricFields. The fieldsCache is per time +// an objectRegistry of all stored fields +template<class GeoField> +static void ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + const word& registryName = "fieldsCache" +); + } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C new file mode 100644 index 0000000000000000000000000000000000000000..30fc2fc0f69fe8fb9018dfc057788f6ffd9013ad --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "interpolationWeights.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +defineTypeNameAndDebug(interpolationWeights, 0); +defineRunTimeSelectionTable(interpolationWeights, word); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +interpolationWeights::interpolationWeights +( + const scalarField& samples +) +: + samples_(samples) +{} + + +// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * // + +autoPtr<interpolationWeights> interpolationWeights::New +( + const word& type, + const scalarField& samples +) +{ + Info<< nl << "Selecting interpolationWeights " + << type << endl; + + wordConstructorTable::iterator cstrIter = + wordConstructorTablePtr_->find(type); + + if (cstrIter == wordConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "interpolationWeights::New(const word&, " + "const scalarField&)" + ) << "Unknown interpolationWeights type " + << type + << endl << endl + << "Valid interpolationWeights types are :" << endl + << wordConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr<interpolationWeights>(cstrIter()(samples)); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +interpolationWeights::~interpolationWeights() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +//objectRegistry& interpolationWeights::registry +//( +// const objectRegistry& obr, +// const word& name +//) +//{ +// if (!obr.foundObject<objectRegistry>(name)) +// { +// objectRegistry* fieldsCachePtr = new objectRegistry +// ( +// IOobject +// ( +// name, +// obr.time().constant(), +// obr, +// IOobject::NO_READ, +// IOobject::NO_WRITE +// ) +// ); +// fieldsCachePtr->store(); +// } +// return const_cast<objectRegistry&>(obr.subRegistry(name)); +//} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H new file mode 100644 index 0000000000000000000000000000000000000000..36e34728356902bd958b4958660cf1e68973f34f --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H @@ -0,0 +1,159 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::interpolationWeights + +Description + Abstract base class for interpolating in 1D + +SourceFiles + interpolationWeights.C + interpolationWeightsTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef interpolationWeights_H +#define interpolationWeights_H + +#include "scalarField.H" +#include "autoPtr.H" +#include "runTimeSelectionTables.H" +#include "pointField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class fvMesh; +class objectRegistry; + +/*---------------------------------------------------------------------------*\ + Class interpolationWeights Declaration +\*---------------------------------------------------------------------------*/ + +class interpolationWeights +{ + +private: + + // Private Member Functions + + //- Disallow default bitwise copy construct + interpolationWeights(const interpolationWeights&); + + //- Disallow default bitwise assignment + void operator=(const interpolationWeights&); + +protected: + + const scalarField& samples_; + +public: + + //- Runtime type information + TypeName("interpolationWeights"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + interpolationWeights, + word, + ( + const scalarField& samples + ), + (samples) + ); + + + // Constructors + + //- Construct from components + interpolationWeights(const scalarField& samples); + + + // Selectors + + //- Return a reference to the selected interpolationWeights + static autoPtr<interpolationWeights> New + ( + const word& type, + const scalarField& samples + ); + + + //- Destructor + virtual ~interpolationWeights(); + + + // Member Functions + + //- Calculate weights and indices to calculate t from samples. + // Returns true if indices changed. + virtual bool valueWeights + ( + const scalar t, + labelList& indices, + scalarField& weights + ) const = 0; + + //- Calculate weights and indices to calculate integrand of t1..t2 + // from samples. Returns true if indices changed. + virtual bool integrationWeights + ( + const scalar t1, + const scalar t2, + labelList& indices, + scalarField& weights + ) const = 0; + + //- Helper: weighted sum + template<class ListType1, class ListType2> + static typename outerProduct + < + typename ListType1::value_type, + typename ListType2::value_type + >::type + weightedSum(const ListType1& f1, const ListType2& f2); + +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "interpolationWeightsTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..9050921b9e59214a2d8cce18205faa645e2f0368 --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C @@ -0,0 +1,77 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "interpolationWeights.H" +#include "ListOps.H" +#include "IOobject.H" +#include "HashSet.H" +#include "objectRegistry.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class ListType1, class ListType2> +typename Foam::outerProduct +< + typename ListType1::value_type, + typename ListType2::value_type +>::type +Foam::interpolationWeights::weightedSum +( + const ListType1& f1, + const ListType2& f2 +) +{ + typedef typename outerProduct + < + typename ListType1::value_type, + typename ListType2::value_type + >::type returnType; + + if (f1.size()) + { + returnType SumProd = f1[0]*f2[0]; + for (label i = 1; i < f1.size(); ++i) + { + SumProd += f1[i]*f2[i]; + } + return SumProd; + } + else + { + return pTraits<returnType>::zero; + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C new file mode 100644 index 0000000000000000000000000000000000000000..18233f1825e64af147de401e28a0e99606cfad11 --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "linearInterpolationWeights.H" +#include "addToRunTimeSelectionTable.H" +#include "ListOps.H" +#include "Pair.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(linearInterpolationWeights, 0); +addToRunTimeSelectionTable +( + interpolationWeights, + linearInterpolationWeights, + word +); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::Pair<Foam::scalar> linearInterpolationWeights::integrationWeights +( + const label i, + const scalar t +) const +{ + // t is in range samples_[i] .. samples_[i+1] + + scalar s = (t-samples_[i])/(samples_[i+1]-samples_[i]); + + if (s < -SMALL || s > 1+SMALL) + { + FatalErrorIn("linearInterpolationWeights::integrationWeights(..)") + << "Value " << t << " outside range " << samples_[i] + << " .. " << samples_[i+1] + << exit(FatalError); + } + + scalar d = samples_[i+1]-t; + + return Pair<scalar>(d*0.5*(1-s), d*0.5*(1+s)); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +linearInterpolationWeights::linearInterpolationWeights +( + const scalarField& samples +) +: + interpolationWeights(samples), + index_(-1) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool linearInterpolationWeights::valueWeights +( + const scalar t, + labelList& indices, + scalarField& weights +) const +{ + bool indexChanged = false; + + // Check if current timeIndex is still valid + if + ( + index_ >= 0 + && index_ < samples_.size() + && ( + samples_[index_] <= t + && (index_ == samples_.size()-1 || t <= samples_[index_+1]) + ) + ) + { + // index_ still at correct slot + } + else + { + // search for correct index + index_ = findLower(samples_, t); + indexChanged = true; + } + + + if (index_ == -1) + { + // Use first element only + indices.setSize(1); + weights.setSize(1); + indices[0] = 0; + weights[0] = 1.0; + } + else if (index_ == samples_.size()-1) + { + // Use last element only + indices.setSize(1); + weights.setSize(1); + indices[0] = samples_.size()-1; + weights[0] = 1.0; + } + else + { + // Interpolate + indices.setSize(2); + weights.setSize(2); + + indices[0] = index_; + indices[1] = index_+1; + + scalar t0 = samples_[indices[0]]; + scalar t1 = samples_[indices[1]]; + scalar deltaT = t1-t0; + + weights[0] = (t1-t)/deltaT; + weights[1] = 1.0-weights[0]; + } + + return indexChanged; +} + + +bool linearInterpolationWeights::integrationWeights +( + const scalar t1, + const scalar t2, + labelList& indices, + scalarField& weights +) const +{ + if (t2 < t1-VSMALL) + { + FatalErrorIn("linearInterpolationWeights::integrationWeights(..)") + << "Integration should be in positive direction." + << " t1:" << t1 << " t2:" << t2 + << exit(FatalError); + } + + // Currently no fancy logic on cached index + label i1 = findLower(samples_, t1); + label i2 = findLower(samples_, t2); + + // For now just fail if any outside table + if (i1 == -1 || i2 == samples_.size()-1) + { + FatalErrorIn("linearInterpolationWeights::integrationWeights(..)") + << "Integrating outside table " << samples_[0] << ".." + << samples_.last() << " not implemented." + << " t1:" << t1 << " t2:" << t2 << exit(FatalError); + } + + label nIndices = i2-i1+2; + + + // Determine if indices already correct + bool anyChanged = false; + + if (nIndices != indices.size()) + { + anyChanged = true; + } + else + { + // Closer check + + label index = i1; + forAll(indices, i) + { + if (indices[i] != index) + { + anyChanged = true; + break; + } + index++; + } + } + + indices.setSize(nIndices); + weights.setSize(nIndices); + weights = 0.0; + + // Sum from i1+1 to i2+1 + for (label i = i1+1; i <= i2; i++) + { + scalar d = samples_[i+1]-samples_[i]; + indices[i-i1] = i; + weights[i-i1] += 0.5*d; + indices[i+1-i1] = i+1; + weights[i+1-i1] += 0.5*d; + } + + // Add from i1 to t1 + { + Pair<scalar> i1Tot1 = integrationWeights(i1, t1); + indices[0] = i1; + weights[0] += i1Tot1.first(); + indices[1] = i1+1; + weights[1] += i1Tot1.second(); + } + + // Subtract from t2 to i2+1 + { + Pair<scalar> wghts = integrationWeights(i2, t2); + indices[i2-i1] = i2; + weights[i2-i1] += -wghts.first(); + indices[i2-i1+1] = i2+1; + weights[i2-i1+1] += -wghts.second(); + } + + return anyChanged; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H new file mode 100644 index 0000000000000000000000000000000000000000..871eb30599cdd5a06d73ad9671af22a190d4152c --- /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 <http://www.gnu.org/licenses/>. + +Class + Foam::linearInterpolationWeights + +Description + +SourceFiles + linearInterpolationWeights.C + +\*---------------------------------------------------------------------------*/ + +#ifndef linearInterpolationWeights_H +#define linearInterpolationWeights_H + +#include "interpolationWeights.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class linearInterpolationWeights Declaration +\*---------------------------------------------------------------------------*/ + +class linearInterpolationWeights +: + public interpolationWeights +{ + +private: + + // Private data + + //- Cached index in samples from previous invocation + mutable label index_; + + // Private Member Functions + + //- Get weights of i and i+1 to calculate integration from t to + // samples_[i+1] + Pair<scalar> integrationWeights + ( + const label i, + const scalar t + ) const; + +public: + + //- Runtime type information + TypeName("linear"); + + // Constructors + + //- Construct from components + linearInterpolationWeights + ( + const scalarField& samples + ); + + + //- Destructor + virtual ~linearInterpolationWeights() + {} + + + // Member Functions + + //- Calculate weights and indices to calculate t from samples. + // Returns true if indices changed. + virtual bool valueWeights + ( + const scalar t, + labelList& indices, + scalarField& weights + ) const; + + //- Calculate weights and indices to calculate integrand of t1..t2 + // from samples. Returns true if indices changed. + virtual bool integrationWeights + ( + const scalar t1, + const scalar t2, + labelList& indices, + scalarField& weights + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C new file mode 100644 index 0000000000000000000000000000000000000000..854146cdff53e41efdcab8204dcdee1883d962e6 --- /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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "splineInterpolationWeights.H" +#include "addToRunTimeSelectionTable.H" +#include "ListOps.H" +#include "linearInterpolationWeights.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(splineInterpolationWeights, 0); +addToRunTimeSelectionTable +( + interpolationWeights, + splineInterpolationWeights, + word +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +splineInterpolationWeights::splineInterpolationWeights +( + const scalarField& samples, + const bool checkEqualDistance +) +: + interpolationWeights(samples), + index_(-1) +{ + if (checkEqualDistance && samples_.size() > 2) + { + const scalar interval = samples_[1]-samples[0]; + for (label i = 2; i < samples_.size(); i++) + { + scalar d = samples_[i]-samples[i-1]; + + if (mag(d-interval) > SMALL) + { + WarningIn + ( + "splineInterpolationWeights::splineInterpolationWeights" + "(const scalarField&)" + ) << "Spline interpolation only valid for constant intervals." + << nl + << "Interval 0-1 : " << interval << nl + << "Interval " << i-1 << '-' << i << " : " + << d << endl; + } + } + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool splineInterpolationWeights::valueWeights +( + const scalar t, + labelList& indices, + scalarField& weights +) const +{ + bool indexChanged = false; + + // linear interpolation + if (samples_.size() <= 2) + { + return linearInterpolationWeights(samples_).valueWeights + ( + t, + indices, + weights + ); + } + + // Check if current timeIndex is still valid + if + ( + index_ >= 0 + && index_ < samples_.size() + && ( + samples_[index_] <= t + && (index_ == samples_.size()-1 || t <= samples_[index_+1]) + ) + ) + { + // index_ still at correct slot + } + else + { + // search for correct index + index_ = findLower(samples_, t); + indexChanged = true; + } + + + // Clamp if outside table + if (index_ == -1) + { + indices.setSize(1); + weights.setSize(1); + + indices[0] = 0; + weights[0] = 1; + return indexChanged; + } + else if (index_ == samples_.size()-1) + { + indices.setSize(1); + weights.setSize(1); + + indices[0] = samples_.size()-1; + weights[0] = 1; + return indexChanged; + } + + + + label lo = index_; + label hi = index_+1; + + // weighting + scalar mu = (t - samples_[lo])/(samples_[hi] - samples_[lo]); + + scalar w0 = 0.5*(mu*(-1+mu*(2-mu))); // coeff of lo-1 + scalar w1 = 0.5*(2+mu*(mu*(-5 + mu*(3)))); // coeff of lo + scalar w2 = 0.5*(mu*(1 + mu*(4 + mu*(-3)))); // coeff of hi + scalar w3 = 0.5*(mu*mu*(-1 + mu)); // coeff of hi+1 + + if (lo > 0) + { + if (hi < samples_.size()-1) + { + // Four points available + indices.setSize(4); + weights.setSize(4); + + indices[0] = lo-1; + indices[1] = lo; + indices[2] = hi; + indices[3] = hi+1; + + weights[0] = w0; + weights[1] = w1; + weights[2] = w2; + weights[3] = w3; + } + else + { + // No y3 available. Extrapolate: y3=3*y2-y1 + indices.setSize(3); + weights.setSize(3); + + indices[0] = lo-1; + indices[1] = lo; + indices[2] = hi; + + weights[0] = w0; + weights[1] = w1 - w3; + weights[2] = w2 + 2*w3; + } + } + else + { + // No y0 available. Extrapolate: y0=2*y1-y2; + if (hi < samples_.size()-1) + { + indices.setSize(3); + weights.setSize(3); + + indices[0] = lo; + indices[1] = hi; + indices[2] = hi+1; + + weights[0] = w1 + 2*w0; + weights[1] = w2 - w0; + weights[2] = w3; + } + else + { + indices.setSize(2); + weights.setSize(2); + + indices[0] = lo; + indices[1] = hi; + + weights[0] = w1 + 2*w0 - w3; + weights[1] = w2 - w0 + 2*w3; + } + } + + return indexChanged; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H new file mode 100644 index 0000000000000000000000000000000000000000..15660ff90cf5233ee310e68148414d274c9bd979 --- /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 <http://www.gnu.org/licenses/>. + +Class + Foam::splineInterpolationWeights + +Description + Catmull-Rom spline interpolation. + +SourceFiles + splineInterpolationWeights.C + +\*---------------------------------------------------------------------------*/ + +#ifndef splineInterpolationWeights_H +#define splineInterpolationWeights_H + +#include "interpolationWeights.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class splineInterpolationWeights Declaration +\*---------------------------------------------------------------------------*/ + +class splineInterpolationWeights +: + public interpolationWeights +{ + +private: + + // Private data + + //- Cached index in samples from previous invocation + mutable label index_; + +public: + + //- Runtime type information + TypeName("spline"); + + // Constructors + + //- Construct from components. By default make sure samples are + // equidistant. + splineInterpolationWeights + ( + const scalarField& samples, + const bool checkEqualDistance = true + ); + + + //- Destructor + virtual ~splineInterpolationWeights() + {} + + + // Member Functions + + //- Calculate weights and indices to calculate t from samples. + // Returns true if indices changed. + virtual bool valueWeights + ( + const scalar t, + labelList& indices, + scalarField& weights + ) const; + + //- Calculate weights and indices to calculate integrand of t1..t2 + // from samples. Returns true if indices changed. + virtual bool integrationWeights + ( + const scalar t1, + const scalar t2, + labelList& indices, + scalarField& weights + ) const + { + notImplemented + ( + "splineInterpolationWeights::integrationWeights(..)" + ); + return false; + } + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C index c977618e2a79ab5c9223e6e17bcd0903eab0c36a..e59c5ec898d3427c074141e7404c7d71dcb9937d 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<class Type> +const Foam::interpolationWeights& Foam::TableBase<Type>::interpolator() const +{ + if (interpolatorPtr_.empty()) + { + // Re-work table into linear list + tableSamples_.setSize(table_.size()); + forAll(table_, i) + { + tableSamples_[i] = table_[i].first(); + } + interpolatorPtr_ = interpolationWeights::New + ( + interpolationScheme_, + tableSamples_ + ); + } + return interpolatorPtr_(); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class Type> @@ -39,6 +62,10 @@ Foam::TableBase<Type>::TableBase(const word& name, const dictionary& dict) dict.lookupOrDefault<word>("outOfBounds", "clamp") ) ), + interpolationScheme_ + ( + dict.lookupOrDefault<word>("interpolationScheme", "linear") + ), table_(), dimensions_(dimless) {} @@ -49,8 +76,11 @@ Foam::TableBase<Type>::TableBase(const TableBase<Type>& tbl) : name_(tbl.name_), boundsHandling_(tbl.boundsHandling_), + interpolationScheme_(tbl.interpolationScheme_), table_(tbl.table_), - dimensions_(tbl.dimensions_) + dimensions_(tbl.dimensions_), + tableSamples_(tbl.tableSamples_), + interpolatorPtr_(tbl.interpolatorPtr_) {} @@ -307,6 +337,9 @@ void Foam::TableBase<Type>::convertTimeBase(const Time& t) scalar value = table_[i].first(); table_[i].first() = t.userTimeToTime(value); } + + tableSamples_.clear(); + interpolatorPtr_.clear(); } @@ -325,88 +358,104 @@ Type Foam::TableBase<Type>::value(const scalar x) const return table_.last().second(); } - // Find i such that x(i) < xDash < x(i+1) - label i = 0; - while ((table_[i+1].first() < xDash) && (i+1 < table_.size())) + // Use interpolator + interpolator().valueWeights(x, currentIndices_, currentWeights_); + + Type t = currentWeights_[0]*table_[currentIndices_[0]].second(); + for (label i = 1; i < currentIndices_.size(); i++) { - i++; + t += currentWeights_[i]*table_[currentIndices_[i]].second(); } - -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() - ); + return t; + + //// Find i such that x(i) < xDash < x(i+1) + //label i = 0; + //while ((table_[i+1].first() < xDash) && (i+1 < table_.size())) + //{ + // i++; + //} + // + //// Linear interpolation to find value + //return Type + //( + // (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first()) + // * (table_[i+1].second() - table_[i].second()) + // + table_[i].second() + //); } template<class Type> Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const { - // Initialise return value - Type sum = pTraits<Type>::zero; - - // 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++; - } + // Use interpolator + interpolator().integrationWeights(x1, x2, currentIndices_, currentWeights_); - // Find next index less than x2 - label id2 = table_.size() - 1; - while ((table_[id2].first() > x2) && (id2 >= 1)) + Type sum = currentWeights_[0]*table_[currentIndices_[0]].second(); + for (label i = 1; i < currentIndices_.size(); i++) { - id2--; + sum += currentWeights_[i]*table_[currentIndices_[i]].second(); } + return sum; - if ((id1 - id2) == 1) - { - // x1 and x2 lie within 1 interval - sum = 0.5*(value(x1) + value(x2))*(x2 - x1); - } - else - { - // x1 and x2 cross multiple intervals - - // Integrate table body - for (label i=id1; i<id2; i++) - { - sum += - (table_[i].second() + table_[i+1].second()) - * (table_[i+1].first() - table_[i].first()); - } - sum *= 0.5; - - // Add table ends (partial segments) - if (id1 > 0) - { - sum += 0.5 - * (value(x1) + table_[id1].second()) - * (table_[id1].first() - x1); - } - if (id2 < table_.size() - 1) - { - sum += 0.5 - * (table_[id2].second() + value(x2)) - * (x2 - table_[id2].first()); - } - } - return sum; + //// Initialise return value + //Type sum = pTraits<Type>::zero; + // + //// Return zero if out of bounds + //if ((x1 > table_.last().first()) || (x2 < table_[0].first())) + //{ + // return sum; + //} + // + //// Find next index greater than x1 + //label id1 = 0; + //while ((table_[id1].first() < x1) && (id1 < table_.size())) + //{ + // id1++; + //} + // + //// Find next index less than x2 + //label id2 = table_.size() - 1; + //while ((table_[id2].first() > x2) && (id2 >= 1)) + //{ + // id2--; + //} + // + //if ((id1 - id2) == 1) + //{ + // // x1 and x2 lie within 1 interval + // sum = 0.5*(value(x1) + value(x2))*(x2 - x1); + //} + //else + //{ + // // x1 and x2 cross multiple intervals + // + // // Integrate table body + // for (label i=id1; i<id2; i++) + // { + // sum += + // (table_[i].second() + table_[i+1].second()) + // * (table_[i+1].first() - table_[i].first()); + // } + // sum *= 0.5; + // + // // Add table ends (partial segments) + // if (id1 > 0) + // { + // sum += 0.5 + // * (value(x1) + table_[id1].second()) + // * (table_[id1].first() - x1); + // } + // if (id2 < table_.size() - 1) + // { + // sum += 0.5 + // * (table_[id2].second() + value(x2)) + // * (x2 - table_[id2].first()); + // } + //} + // + //return sum; } diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H index 15164b3d903819ebd1a79f25980548238e02e039..ce4c3d2287960883c8f2cc63457d871467dacd0b 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<Tuple2<scalar, Type> > table_; @@ -91,9 +95,24 @@ protected: //- The dimension set dimensionSet dimensions_; + //- Extracted values + mutable scalarField tableSamples_; + + //- Interpolator method + mutable autoPtr<interpolationWeights> interpolatorPtr_; + + //- Cached indices and weights + mutable labelList currentIndices_; + + mutable scalarField currentWeights_; + // Protected Member Functions + + //- Return (demand driven) interpolator + const interpolationWeights& interpolator() const; + //- Disallow default bitwise assignment void operator=(const TableBase<Type>&); diff --git a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C index 5323ac6a859f4590094fe249825b6ec3e63df771..15fe8b57a265ef270198339197f9656200261ee6 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<Type>::TableFile(const word& entryName, const dictionary& dict) fileName expandedFile(fName_); IFstream is(expandedFile.expand()); + if (!is.good()) + { + FatalIOErrorIn + ( + "TableFile<Type>::TableFile(const word&, const dictionary&)", + is + ) << "Cannot open file." << exit(FatalIOError); + } + is >> this->table_; TableBase<Type>::check(); diff --git a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H index f167c211561dc64c0a627a3fee9db7701be1869e..7948ac6ee666be3cb506a1adb2944a2a2c8d28d2 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H @@ -31,9 +31,10 @@ Description <entryName> tableFile; tableFileCoeffs { - dimensions [0 0 1 0 0]; // optional dimensions - fileName dataFile; // name of data file - outOfBounds clamp; // optional out-of-bounds handling + dimensions [0 0 1 0 0]; // optional dimensions + fileName dataFile; // name of data file + outOfBounds clamp; // optional out-of-bounds handling + interpolationScheme linear; // optional interpolation method } \endverbatim diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C index 16704d4adc332c582a3b06557951a19935a2d91e..2d31d3802f1630eef14cb1310a29fde1e80c39f0 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C @@ -209,7 +209,7 @@ template<class Type> void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable() { // Initialise - if (startSampleTime_ == -1 && endSampleTime_ == -1) + if (mapperPtr_.empty()) { pointIOField samplePoints ( diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files index 43dc7513e7dd322d16322e1f6bc82d74d97f9577..cdaf67e6afa5fe670d13d35b24faf34c2c619560 100644 --- a/src/fvMotionSolver/Make/files +++ b/src/fvMotionSolver/Make/files @@ -38,6 +38,7 @@ $(derivedPoint)/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C $(derivedPoint)/waveDisplacement/waveDisplacementPointPatchVectorField.C $(derivedPoint)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C +$(derivedPoint)/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C diff --git a/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000000000000000000000000000000000..21b5e74224fbf6fa79330dbbe9372195f92f9d9c --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C @@ -0,0 +1,291 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "uniformInterpolatedDisplacementPointPatchVectorField.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "polyMesh.H" +#include "interpolationWeights.H" +#include "uniformInterpolate.H" +#include "ReadFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +uniformInterpolatedDisplacementPointPatchVectorField:: +uniformInterpolatedDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF +) +: + fixedValuePointPatchField<vector>(p, iF) +{} + + +uniformInterpolatedDisplacementPointPatchVectorField:: +uniformInterpolatedDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF, + const dictionary& dict +) +: + fixedValuePointPatchField<vector>(p, iF, dict), + fieldName_(dict.lookup("fieldName")), + interpolationScheme_(dict.lookup("interpolationScheme")) +{ + const pointMesh& pMesh = this->dimensionedInternalField().mesh(); + + // Read time values + instantList allTimes = Time::findTimes(pMesh().time().path()); + + // Only keep those that contain the field + DynamicList<word> names(allTimes.size()); + DynamicList<scalar> values(allTimes.size()); + + forAll(allTimes, i) + { + IOobject io + ( + fieldName_, + allTimes[i].name(), + pMesh(), + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ); + if (io.headerOk()) + { + names.append(allTimes[i].name()); + values.append(allTimes[i].value()); + } + } + timeNames_.transfer(names); + timeVals_.transfer(values); + + Info<< type() << " : found " << fieldName_ << " for times " + << timeNames_ << endl; + + if (timeNames_.size() < 1) + { + FatalErrorIn + ( + "uniformInterpolatedDisplacementPointPatchVectorField::\n" + "uniformInterpolatedDisplacementPointPatchVectorField\n" + "(\n" + " const pointPatch&,\n" + " const DimensionedField<vector, pointMesh>&,\n" + " const dictionary&\n" + ")\n" + ) << "Did not find any times with " << fieldName_ + << exit(FatalError); + } + + + if (!dict.found("value")) + { + updateCoeffs(); + } +} + + +uniformInterpolatedDisplacementPointPatchVectorField:: +uniformInterpolatedDisplacementPointPatchVectorField +( + const uniformInterpolatedDisplacementPointPatchVectorField& ptf, + const pointPatch& p, + const DimensionedField<vector, pointMesh>& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField<vector>(ptf, p, iF, mapper), + fieldName_(ptf.fieldName_), + interpolationScheme_(ptf.interpolationScheme_), + timeNames_(ptf.timeNames_), + timeVals_(ptf.timeVals_), + interpolatorPtr_(ptf.interpolatorPtr_) +{} + + +uniformInterpolatedDisplacementPointPatchVectorField:: +uniformInterpolatedDisplacementPointPatchVectorField +( + const uniformInterpolatedDisplacementPointPatchVectorField& ptf, + const DimensionedField<vector, pointMesh>& iF +) +: + fixedValuePointPatchField<vector>(ptf, iF), + fieldName_(ptf.fieldName_), + interpolationScheme_(ptf.interpolationScheme_), + timeNames_(ptf.timeNames_), + timeVals_(ptf.timeVals_), + interpolatorPtr_(ptf.interpolatorPtr_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void uniformInterpolatedDisplacementPointPatchVectorField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + if (!interpolatorPtr_.valid()) + { + interpolatorPtr_ = interpolationWeights::New + ( + interpolationScheme_, + timeVals_ + ); + } + + const pointMesh& pMesh = this->dimensionedInternalField().mesh(); + const Time& t = pMesh().time(); + + // Update indices of times and weights + bool timesChanged = interpolatorPtr_->valueWeights + ( + t.timeOutputValue(), + currentIndices_, + currentWeights_ + ); + + const wordList currentTimeNames + ( + UIndirectList<word>(timeNames_, currentIndices_) + ); + + + // Load if necessary fields for this interpolation + if (timesChanged) + { + objectRegistry& fieldsCache = const_cast<objectRegistry&> + ( + pMesh.thisDb().subRegistry("fieldsCache", true) + ); + // Save old times so we now which ones have been loaded and need + // 'correctBoundaryConditions'. Bit messy. + HashSet<word> oldTimes(fieldsCache.toc()); + + ReadFields<pointVectorField> + ( + fieldName_, + pMesh, + currentTimeNames + ); + + forAllConstIter(objectRegistry, fieldsCache, fieldsCacheIter) + { + if (!oldTimes.found(fieldsCacheIter.key())) + { + // Newly loaded fields. Make sure the internal + // values are consistent with the boundary conditions. + // This is quite often not the case since these + // fields typically are constructed 'by hand' + + const objectRegistry& timeCache = dynamic_cast + < + const objectRegistry& + >(*fieldsCacheIter()); + + + pointVectorField& d = const_cast<pointVectorField&> + ( + timeCache.lookupObject<pointVectorField> + ( + fieldName_ + ) + ); + d.correctBoundaryConditions(); + } + } + } + + + // Interpolate the whole field + pointVectorField result + ( + uniformInterpolate<pointVectorField> + ( + IOobject + ( + word("uniformInterpolate(") + + this->dimensionedInternalField().name() + + ')', + pMesh.time().timeName(), + pMesh.thisDb(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fieldName_, + currentTimeNames, + currentWeights_ + ) + ); + + + // Extract back from the internal field + this->operator== + ( + this->patchInternalField(result.dimensionedInternalField()) + ); + + fixedValuePointPatchField<vector>::updateCoeffs(); +} + + +void uniformInterpolatedDisplacementPointPatchVectorField::write(Ostream& os) +const +{ + pointPatchField<vector>::write(os); + os.writeKeyword("fieldName") + << fieldName_ << token::END_STATEMENT << nl; + os.writeKeyword("interpolationScheme") + << interpolationScheme_ << token::END_STATEMENT << nl; + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchTypeField +( + pointPatchVectorField, + uniformInterpolatedDisplacementPointPatchVectorField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.H b/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.H new file mode 100644 index 0000000000000000000000000000000000000000..bd1beb36ff6523c2dff32677a679c8938f4e8f21 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.H @@ -0,0 +1,183 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::uniformInterpolatedDisplacementPointPatchVectorField + +Description + Interpolates pre-specified motion. + + Motion specified as pointVectorFields. E.g. + + walls + { + type uniformInterpolatedDisplacement; + value uniform (0 0 0); + fieldName wantedDisplacement; + interpolationScheme linear; + } + + This will scan the case for 'wantedDisplacement' pointVectorFields + and interpolate those in time (using 'linear' interpolation) to + obtain the current displacement. + The advantage of specifying displacement in this way is that it + automatically works through decomposePar. + +SourceFiles + uniformInterpolatedDisplacementPointPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef uniformInterpolatedDisplacementPointPatchVectorField_H +#define uniformInterpolatedDisplacementPointPatchVectorField_H + +#include "fixedValuePointPatchField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class interpolationWeights; + +/*---------------------------------------------------------------------------*\ + Class uniformInterpolatedDisplacementPointPatchVectorField Declaration +\*---------------------------------------------------------------------------*/ + +class uniformInterpolatedDisplacementPointPatchVectorField +: + public fixedValuePointPatchField<vector> +{ + // Private data + + //- Name of displacement field + const word fieldName_; + + const word interpolationScheme_; + + //- Times with pre-specified displacement + wordList timeNames_; + + //- Times with pre-specified displacement + scalarField timeVals_; + + //- User-specified interpolator + autoPtr<interpolationWeights> interpolatorPtr_; + + + //- Cached interpolation times + labelList currentIndices_; + + //- Cached interpolation weights + scalarField currentWeights_; + +public: + + //- Runtime type information + TypeName("uniformInterpolatedDisplacement"); + + + // Constructors + + //- Construct from patch and internal field + uniformInterpolatedDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField<vector, pointMesh>& + ); + + //- Construct from patch, internal field and dictionary + uniformInterpolatedDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField<vector, pointMesh>&, + const dictionary& + ); + + //- Construct by mapping given patchField<vector> onto a new patch + uniformInterpolatedDisplacementPointPatchVectorField + ( + const uniformInterpolatedDisplacementPointPatchVectorField&, + const pointPatch&, + const DimensionedField<vector, pointMesh>&, + const pointPatchFieldMapper& + ); + + //- Construct and return a clone + virtual autoPtr<pointPatchField<vector> > clone() const + { + return autoPtr<pointPatchField<vector> > + ( + new uniformInterpolatedDisplacementPointPatchVectorField + ( + *this + ) + ); + } + + //- Construct as copy setting internal field reference + uniformInterpolatedDisplacementPointPatchVectorField + ( + const uniformInterpolatedDisplacementPointPatchVectorField&, + const DimensionedField<vector, pointMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr<pointPatchField<vector> > clone + ( + const DimensionedField<vector, pointMesh>& iF + ) const + { + return autoPtr<pointPatchField<vector> > + ( + new uniformInterpolatedDisplacementPointPatchVectorField + ( + *this, + iF + ) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C index ccf2b0c38254eba351d3a3f10b5113f591a151f2..aedaec8bcea7509f5aecac45c93c235a0c3aa5ef 100644 --- a/src/surfMesh/MeshedSurface/MeshedSurface.C +++ b/src/surfMesh/MeshedSurface/MeshedSurface.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 @@ -656,6 +656,7 @@ bool Foam::MeshedSurface<Face>::stitchFaces << " faces" << endl; } faceLst.setSize(newFaceI); + faceMap.setSize(newFaceI); remapFaces(faceMap); } faceMap.clear();