diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C index 61fb24ba77d505f1bd996c4bec52f4563978a259..1831bd903ddfc5d1a5e98db9432aa19f01d749ca 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.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 @@ -25,12 +25,7 @@ License #include "timeVaryingMappedFixedValueFvPatchField.H" #include "Time.H" -#include "triSurfaceTools.H" -#include "triSurface.H" -#include "vector2D.H" -#include "OFstream.H" #include "AverageIOField.H" -#include "Random.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -51,9 +46,6 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(iF.name()), setAverage_(false), perturb_(0), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -78,9 +70,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -104,9 +94,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(iF.name()), setAverage_(readBool(dict.lookup("setAverage"))), perturb_(dict.lookupOrDefault("perturb", 1E-5)), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), + mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -139,9 +127,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(ptf.referenceCS_), - nearestVertex_(ptf.nearestVertex_), - nearestVertexWeight_(ptf.nearestVertexWeight_), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), startSampledValues_(ptf.startSampledValues_), @@ -165,9 +151,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(ptf.referenceCS_), - nearestVertex_(ptf.nearestVertex_), - nearestVertexWeight_(ptf.nearestVertexWeight_), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), startSampledValues_(ptf.startSampledValues_), @@ -192,6 +176,10 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::autoMap startSampledValues_.autoMap(m); endSampledValues_.autoMap(m); } + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; } @@ -209,266 +197,90 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::rmap startSampledValues_.rmap(tiptf.startSampledValues_, addr); endSampledValues_.rmap(tiptf.endSampledValues_, addr); + + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; } template<class Type> -void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints() +void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable() { - // Read the sample points - - pointIOField samplePoints - ( - IOobject - ( - "points", - this->db().time().constant(), - "boundaryData"/this->patch().name(), - this->db(), - IOobject::MUST_READ, - IOobject::AUTO_WRITE, - false - ) - ); - - const fileName samplePointsFile = samplePoints.filePath(); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Read " << samplePoints.size() << " sample points from " - << samplePointsFile << endl; - } - - // Determine coordinate system from samplePoints - - if (samplePoints.size() < 3) + // Initialise + if (startSampleTime_ == -1 && endSampleTime_ == -1) { - FatalErrorIn + pointIOField samplePoints ( - "timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()" - ) << "Only " << samplePoints.size() << " points read from file " - << samplePoints.objectPath() << nl - << "Need at least three non-colinear samplePoints" - << " to be able to interpolate." - << "\n on patch " << this->patch().name() - << " of points " << samplePoints.name() - << " in file " << samplePoints.objectPath() - << exit(FatalError); - } - - const point& p0 = samplePoints[0]; - - // Find furthest away point - vector e1; - label index1 = -1; - scalar maxDist = -GREAT; - - for (label i = 1; i < samplePoints.size(); i++) - { - const vector d = samplePoints[i] - p0; - scalar magD = mag(d); + IOobject + ( + "points", + this->db().time().constant(), + "boundaryData"/this->patch().name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + const fileName samplePointsFile = samplePoints.filePath(); - if (magD > maxDist) + if (debug) { - e1 = d/magD; - index1 = i; - maxDist = magD; + Info<< "timeVaryingMappedFixedValueFvPatchField :" + << " Read " << samplePoints.size() << " sample points from " + << samplePointsFile << endl; } - } - // Find point that is furthest away from line p0-p1 - const point& p1 = samplePoints[index1]; - label index2 = -1; - maxDist = -GREAT; - for (label i = 1; i < samplePoints.size(); i++) - { - if (i != index1) - { - const point& p2 = samplePoints[i]; - vector e2(p2 - p0); - e2 -= (e2&e1)*e1; - scalar magE2 = mag(e2); - - if (magE2 > maxDist) - { - index2 = i; - maxDist = magE2; - } - } - } - if (index2 == -1) - { - FatalErrorIn + // Allocate the interpolator + mapperPtr_.reset ( - "timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()" - ) << "Cannot find points that make valid normal." << nl - << "Have so far points " << p0 << " and " << p1 - << "Need at least three sample points which are not in a line." - << "\n on patch " << this->patch().name() - << " of points " << samplePoints.name() - << " in file " << samplePoints.objectPath() - << exit(FatalError); - } - - vector n = e1^(samplePoints[index2]-p0); - n /= mag(n); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Used points " << p0 << ' ' << samplePoints[index1] - << ' ' << samplePoints[index2] - << " to define coordinate system with normal " << n << endl; - } - - referenceCS_.reset - ( - new coordinateSystem - ( - "reference", - p0, // origin - n, // normal - e1 // 0-axis - ) - ); - - tmp<vectorField> tlocalVertices - ( - referenceCS().localPosition(samplePoints) - ); - vectorField& localVertices = tlocalVertices(); - - const boundBox bb(localVertices, true); - const point bbMid(bb.midpoint()); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Perturbing points with " << perturb_ - << " fraction of a random position inside " << bb - << " to break any ties on regular meshes." - << nl << endl; - } - - Random rndGen(123456); - forAll(localVertices, i) - { - localVertices[i] += - perturb_ - *(rndGen.position(bb.min(), bb.max())-bbMid); - } - - // Determine triangulation - List<vector2D> localVertices2D(localVertices.size()); - forAll(localVertices, i) - { - localVertices2D[i][0] = localVertices[i][0]; - localVertices2D[i][1] = localVertices[i][1]; - } - - triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); - - tmp<pointField> tlocalFaceCentres - ( - referenceCS().localPosition - ( - this->patch().patch().faceCentres() - ) - ); - const pointField& localFaceCentres = tlocalFaceCentres(); - - if (debug) - { - Pout<< "readSamplePoints :" - <<" Dumping triangulated surface to triangulation.stl" << endl; - s.write(this->db().time().path()/"triangulation.stl"); + new pointToPointPlanarInterpolation + ( + samplePoints, + this->patch().patch().faceCentres(), + perturb_ + ) + ); - OFstream str(this->db().time().path()/"localFaceCentres.obj"); - Pout<< "readSamplePoints :" - << " Dumping face centres to " << str.name() << endl; + // Read the times for which data is available + const fileName samplePointsDir = samplePointsFile.path(); + sampleTimes_ = Time::findTimes(samplePointsDir); - forAll(localFaceCentres, i) + if (debug) { - const point& p = localFaceCentres[i]; - str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; + Info<< "timeVaryingMappedFixedValueFvPatchField : In directory " + << samplePointsDir << " found times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) + << endl; } } - // Determine interpolation onto face centres. - triSurfaceTools::calcInterpolationWeights - ( - s, - localFaceCentres, // points to interpolate to - nearestVertex_, - nearestVertexWeight_ - ); - - - - // Read the times for which data is available - - const fileName samplePointsDir = samplePointsFile.path(); - - sampleTimes_ = Time::findTimes(samplePointsDir); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField : In directory " - << samplePointsDir << " found times " << timeNames(sampleTimes_) - << endl; - } -} - - -template<class Type> -wordList timeVaryingMappedFixedValueFvPatchField<Type>::timeNames -( - const instantList& times -) -{ - wordList names(times.size()); - - forAll(times, i) - { - names[i] = times[i].name(); - } - return names; -} - -template<class Type> -void timeVaryingMappedFixedValueFvPatchField<Type>::findTime -( - const fileName& instance, - const fileName& local, - const scalar timeVal, - label& lo, - label& hi -) const -{ - lo = startSampleTime_; - hi = -1; + // Find current time in sampleTimes + label lo = -1; + label hi = -1; - for (label i = startSampleTime_+1; i < sampleTimes_.size(); i++) - { - if (sampleTimes_[i].value() > timeVal) - { - break; - } - else - { - lo = i; - } - } + bool foundTime = mapperPtr_().findTime + ( + sampleTimes_, + startSampleTime_, + this->db().time().value(), + lo, + hi + ); - if (lo == -1) + if (!foundTime) { - FatalErrorIn("findTime") - << "Cannot find starting sampling values for current time " - << timeVal << nl + FatalErrorIn + ( + "timeVaryingMappedFixedValueFvPatchField<Type>::checkTable" + ) << "Cannot find starting sampling values for current time " + << this->db().time().value() << nl << "Have sampling values for times " - << timeNames(sampleTimes_) << nl + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl << "In directory " << this->db().time().constant()/"boundaryData"/this->patch().name() << "\n on patch " << this->patch().name() @@ -476,52 +288,6 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::findTime << exit(FatalError); } - if (lo < sampleTimes_.size()-1) - { - hi = lo+1; - } - - - if (debug) - { - if (hi == -1) - { - Pout<< "findTime : Found time " << timeVal << " after" - << " index:" << lo << " time:" << sampleTimes_[lo].value() - << endl; - } - else - { - Pout<< "findTime : Found time " << timeVal << " inbetween" - << " index:" << lo << " time:" << sampleTimes_[lo].value() - << " and index:" << hi << " time:" << sampleTimes_[hi].value() - << endl; - } - } -} - - -template<class Type> -void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable() -{ - // Initialise - if (startSampleTime_ == -1 && endSampleTime_ == -1) - { - readSamplePoints(); - } - - // Find current time in sampleTimes - label lo = -1; - label hi = -1; - - findTime - ( - this->db().time().constant(), - "boundaryData"/this->patch().name(), - this->db().time().value(), - lo, - hi - ); // Update sampled data fields. @@ -573,7 +339,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable() ); startAverage_ = vals.average(); - startSampledValues_ = interpolate(vals); + startSampledValues_ = mapperPtr_().interpolate(vals); } } @@ -617,61 +383,22 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable() ) ); endAverage_ = vals.average(); - endSampledValues_ = interpolate(vals); + endSampledValues_ = mapperPtr_().interpolate(vals); } } } -template<class Type> -tmp<Field<Type> > timeVaryingMappedFixedValueFvPatchField<Type>::interpolate -( - const Field<Type>& sourceFld -) const -{ - tmp<Field<Type> > tfld(new Field<Type>(nearestVertex_.size())); - Field<Type>& fld = tfld(); - - forAll(fld, i) - { - const FixedList<label, 3>& verts = nearestVertex_[i]; - const FixedList<scalar, 3>& w = nearestVertexWeight_[i]; - - if (verts[2] == -1) - { - if (verts[1] == -1) - { - // Use vertex0 only - fld[i] = sourceFld[verts[0]]; - } - else - { - // Use vertex 0,1 - fld[i] = - w[0]*sourceFld[verts[0]] - + w[1]*sourceFld[verts[1]]; - } - } - else - { - fld[i] = - w[0]*sourceFld[verts[0]] - + w[1]*sourceFld[verts[1]] - + w[2]*sourceFld[verts[2]]; - } - } - return tfld; -} - - template<class Type> void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs() { + if (this->updated()) { return; } + checkTable(); // Interpolate between the sampled data @@ -768,7 +495,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const { fvPatchField<Type>::write(os); os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl; - os.writeKeyword("peturb") << perturb_ << token::END_STATEMENT << nl; + os.writeKeyword("perturb") << perturb_ << token::END_STATEMENT << nl; if (fieldTableName_ != this->dimensionedInternalField().name()) { diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H index e3ac11215495207ba2572e0a8f978434c7a913f8..1494efec357831f721c950156d06a7c0795d199f 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.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 @@ -66,9 +66,9 @@ SourceFiles #define timeVaryingMappedFixedValueFvPatchField_H #include "fixedValueFvPatchFields.H" -#include "coordinateSystem.H" #include "FixedList.H" #include "instantList.H" +#include "pointToPointPlanarInterpolation.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -95,16 +95,8 @@ class timeVaryingMappedFixedValueFvPatchField //- Fraction of perturbation (fraction of bounding box) to add scalar perturb_; - //- Coordinate system - autoPtr<coordinateSystem> referenceCS_; - - //- Current interpolation addressing to face centres of underlying - // patch - List<FixedList<label, 3> > nearestVertex_; - - //- Current interpolation factors to face centres of underlying - // patch - List<FixedList<scalar, 3> > nearestVertexWeight_; + //- 2D interpolation + autoPtr<pointToPointPlanarInterpolation> mapperPtr_; //- List of boundaryData time directories instantList sampleTimes_; @@ -127,31 +119,6 @@ class timeVaryingMappedFixedValueFvPatchField //- If setAverage: end average value Type endAverage_; - - // Private Member Functions - - //- Get names of times - static wordList timeNames(const instantList&); - - //- Find times around current time - void findTime - ( - const fileName& instance, - const fileName& local, - const scalar timeVal, - label& lo, - label& hi - ) const; - - - //- Read boundary points and determine interpolation weights to patch - // faceCentres - void readSamplePoints(); - - //- Do actual interpolation using current weights - tmp<Field<Type> > interpolate(const Field<Type>&) const; - - public: //- Runtime type information @@ -224,12 +191,6 @@ public: // Access - //- Return the coordinateSystem - const coordinateSystem& referenceCS() const - { - return referenceCS_; - } - //- Return startSampledValues const Field<Type> startSampledValues() { diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files index bea2f1065699047d4470169b4a9b92cbb1e1e0c7..43dc7513e7dd322d16322e1f6bc82d74d97f9577 100644 --- a/src/fvMotionSolver/Make/files +++ b/src/fvMotionSolver/Make/files @@ -26,13 +26,19 @@ motionDiffusivity/manipulators/exponential/exponentialDiffusivity.C fvPatchFields/derived/cellMotion/cellMotionFvPatchFields.C fvPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementFvPatchFields.C -pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C -pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C - -pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C -pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C -pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C -pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C -pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C +derivedPoint = pointPatchFields/derived + +$(derivedPoint)/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C +$(derivedPoint)/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C + +$(derivedPoint)/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C +$(derivedPoint)/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C +$(derivedPoint)/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C +$(derivedPoint)/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C +$(derivedPoint)/waveDisplacement/waveDisplacementPointPatchVectorField.C + +$(derivedPoint)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C + + LIB = $(FOAM_LIBBIN)/libfvMotionSolvers diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C new file mode 100644 index 0000000000000000000000000000000000000000..e1bdfc0910262ec9f9d24adf1d06a9e577b4988b --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C @@ -0,0 +1,482 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "timeVaryingMappedFixedValuePointPatchField.H" +#include "Time.H" +#include "AverageIOField.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class Type> +Foam:: +timeVaryingMappedFixedValuePointPatchField<Type>:: +timeVaryingMappedFixedValuePointPatchField +( + const pointPatch& p, + const DimensionedField<Type, pointMesh>& iF +) +: + fixedValuePointPatchField<Type>(p, iF), + fieldTableName_(iF.name()), + setAverage_(false), + perturb_(0), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits<Type>::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits<Type>::zero) + +{} + + +template<class Type> +Foam:: +timeVaryingMappedFixedValuePointPatchField<Type>:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField<Type>& ptf, + const pointPatch& p, + const DimensionedField<Type, pointMesh>& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField<Type>(ptf, p, iF, mapper), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits<Type>::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits<Type>::zero) +{} + + +template<class Type> +Foam:: +timeVaryingMappedFixedValuePointPatchField<Type>:: +timeVaryingMappedFixedValuePointPatchField +( + const pointPatch& p, + const DimensionedField<Type, pointMesh>& iF, + const dictionary& dict +) +: + fixedValuePointPatchField<Type>(p, iF), + fieldTableName_(iF.name()), + setAverage_(readBool(dict.lookup("setAverage"))), + perturb_(dict.lookupOrDefault("perturb", 1E-5)), + mapperPtr_(NULL), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits<Type>::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits<Type>::zero) +{ + dict.readIfPresent("fieldTableName", fieldTableName_); + updateCoeffs(); +} + + +template<class Type> +Foam:: +timeVaryingMappedFixedValuePointPatchField<Type>:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField<Type>& ptf +) +: + fixedValuePointPatchField<Type>(ptf), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(ptf.sampleTimes_), + startSampleTime_(ptf.startSampleTime_), + startSampledValues_(ptf.startSampledValues_), + startAverage_(ptf.startAverage_), + endSampleTime_(ptf.endSampleTime_), + endSampledValues_(ptf.endSampledValues_), + endAverage_(ptf.endAverage_) +{} + + +template<class Type> +Foam:: +timeVaryingMappedFixedValuePointPatchField<Type>:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField<Type>& ptf, + const DimensionedField<Type, pointMesh>& iF +) +: + fixedValuePointPatchField<Type>(ptf, iF), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(ptf.sampleTimes_), + startSampleTime_(ptf.startSampleTime_), + startSampledValues_(ptf.startSampledValues_), + startAverage_(ptf.startAverage_), + endSampleTime_(ptf.endSampleTime_), + endSampledValues_(ptf.endSampledValues_), + endAverage_(ptf.endAverage_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable() +{ + // Initialise + if (startSampleTime_ == -1 && endSampleTime_ == -1) + { + const polyMesh& pMesh = this->patch().boundaryMesh().mesh()(); + + // Read the initial point position + pointField meshPts; + + if (pMesh.pointsInstance() == pMesh.facesInstance()) + { + meshPts = pointField(pMesh.points(), this->patch().meshPoints()); + } + else + { + // Load points from facesInstance + if (debug) + { + Info<< "Reloading points0 from " << pMesh.facesInstance() + << endl; + } + + pointIOField points0 + ( + IOobject + ( + "points", + pMesh.facesInstance(), + polyMesh::meshSubDir, + pMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) + ); + meshPts = pointField(points0, this->patch().meshPoints()); + } + + pointIOField samplePoints + ( + IOobject + ( + "points", + this->db().time().constant(), + "boundaryData"/this->patch().name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + mapperPtr_.reset + ( + new pointToPointPlanarInterpolation + ( + samplePoints, + meshPts, + perturb_ + ) + ); + + // Read the times for which data is available + + const fileName samplePointsFile = samplePoints.filePath(); + const fileName samplePointsDir = samplePointsFile.path(); + sampleTimes_ = Time::findTimes(samplePointsDir); + + if (debug) + { + Info<< "timeVaryingMappedFixedValuePointPatchField : In directory " + << samplePointsDir << " found times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) + << endl; + } + } + + // Find current time in sampleTimes + label lo = -1; + label hi = -1; + + bool foundTime = mapperPtr_().findTime + ( + sampleTimes_, + startSampleTime_, + this->db().time().value(), + lo, + hi + ); + + if (!foundTime) + { + FatalErrorIn + ( + "timeVaryingMappedFixedValuePointPatchField<Type>::checkTable" + ) << "Cannot find starting sampling values for current time " + << this->db().time().value() << nl + << "Have sampling values for times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl + << "In directory " + << this->db().time().constant()/"boundaryData"/this->patch().name() + << "\n on patch " << this->patch().name() + << " of field " << fieldTableName_ + << exit(FatalError); + } + + + // Update sampled data fields. + + if (lo != startSampleTime_) + { + startSampleTime_ = lo; + + if (startSampleTime_ == endSampleTime_) + { + // No need to reread since are end values + if (debug) + { + Pout<< "checkTable : Setting startValues to (already read) " + << "boundaryData" + /this->patch().name() + /sampleTimes_[startSampleTime_].name() + << endl; + } + startSampledValues_ = endSampledValues_; + startAverage_ = endAverage_; + } + else + { + if (debug) + { + Pout<< "checkTable : Reading startValues from " + << "boundaryData" + /this->patch().name() + /sampleTimes_[lo].name() + << endl; + } + + // Reread values and interpolate + AverageIOField<Type> vals + ( + IOobject + ( + fieldTableName_, + this->db().time().constant(), + "boundaryData" + /this->patch().name() + /sampleTimes_[startSampleTime_].name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + startAverage_ = vals.average(); + startSampledValues_ = mapperPtr_().interpolate(vals); + } + } + + if (hi != endSampleTime_) + { + endSampleTime_ = hi; + + if (endSampleTime_ == -1) + { + // endTime no longer valid. Might as well clear endValues. + if (debug) + { + Pout<< "checkTable : Clearing endValues" << endl; + } + endSampledValues_.clear(); + } + else + { + if (debug) + { + Pout<< "checkTable : Reading endValues from " + << "boundaryData" + /this->patch().name() + /sampleTimes_[endSampleTime_].name() + << endl; + } + // Reread values and interpolate + AverageIOField<Type> vals + ( + IOobject + ( + fieldTableName_, + this->db().time().constant(), + "boundaryData" + /this->patch().name() + /sampleTimes_[endSampleTime_].name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + endAverage_ = vals.average(); + endSampledValues_ = mapperPtr_().interpolate(vals); + } + } +} + + +template<class Type> +void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + checkTable(); + + // Interpolate between the sampled data + + Type wantedAverage; + + if (endSampleTime_ == -1) + { + // only start value + if (debug) + { + Pout<< "updateCoeffs : Sampled, non-interpolated values" + << " from start time:" + << sampleTimes_[startSampleTime_].name() << nl; + } + + this->operator==(startSampledValues_); + wantedAverage = startAverage_; + } + else + { + scalar start = sampleTimes_[startSampleTime_].value(); + scalar end = sampleTimes_[endSampleTime_].value(); + + scalar s = (this->db().time().value()-start)/(end-start); + + if (debug) + { + Pout<< "updateCoeffs : Sampled, interpolated values" + << " between start time:" + << sampleTimes_[startSampleTime_].name() + << " and end time:" << sampleTimes_[endSampleTime_].name() + << " with weight:" << s << endl; + } + + this->operator==((1-s)*startSampledValues_ + s*endSampledValues_); + wantedAverage = (1-s)*startAverage_ + s*endAverage_; + } + + // Enforce average. Either by scaling (if scaling factor > 0.5) or by + // offsetting. + if (setAverage_) + { + const Field<Type>& fld = *this; + + Type averagePsi = gAverage(fld); + + if (debug) + { + Pout<< "updateCoeffs :" + << " actual average:" << averagePsi + << " wanted average:" << wantedAverage + << endl; + } + + if (mag(averagePsi) < VSMALL) + { + // Field too small to scale. Offset instead. + const Type offset = wantedAverage - averagePsi; + if (debug) + { + Pout<< "updateCoeffs :" + << " offsetting with:" << offset << endl; + } + this->operator==(fld+offset); + } + else + { + const scalar scale = mag(wantedAverage)/mag(averagePsi); + + if (debug) + { + Pout<< "updateCoeffs :" + << " scaling with:" << scale << endl; + } + this->operator==(scale*fld); + } + } + + if (debug) + { + Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this) + << " max:" << gMax(*this) << endl; + } + + fixedValuePointPatchField<Type>::updateCoeffs(); +} + + +template<class Type> +void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::write +( + Ostream& os +) const +{ + fixedValuePointPatchField<Type>::write(os); + os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl; + os.writeKeyword("perturb") << perturb_ << token::END_STATEMENT << nl; + + if (fieldTableName_ != this->dimensionedInternalField().name()) + { + os.writeKeyword("fieldTableName") << fieldTableName_ + << token::END_STATEMENT << nl; + } +} + + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H new file mode 100644 index 0000000000000000000000000000000000000000..108ec6786ee9b4a804ec0e5c1ebe0068f56ae538 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H @@ -0,0 +1,195 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::timeVaryingMappedFixedValuePointPatchField + +Description + A time-varying form of a mapped fixed value boundary condition. + +See Also + Foam::timeVaryingMappedFixedValueFvPatchField + +SourceFiles + timeVaryingMappedFixedValuePointPatchField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingMappedFixedValuePointPatchField_H +#define timeVaryingMappedFixedValuePointPatchField_H + +#include "fixedValuePointPatchField.H" +#include "instantList.H" +#include "pointToPointPlanarInterpolation.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class timeVaryingMappedFixedValuePointPatchField Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type> +class timeVaryingMappedFixedValuePointPatchField +: + public fixedValuePointPatchField<Type> +{ + // Private data + + //- Name of the field data table, defaults to the name of the field + word fieldTableName_; + + //- If true adjust the mapped field to maintain average value + bool setAverage_; + + //- Fraction of perturbation (fraction of bounding box) to add + scalar perturb_; + + //- 2D interpolation + autoPtr<pointToPointPlanarInterpolation> mapperPtr_; + + //- List of boundaryData time directories + instantList sampleTimes_; + + //- Current starting index in sampleTimes + label startSampleTime_; + + //- Interpolated values from startSampleTime + Field<Type> startSampledValues_; + + //- If setAverage: starting average value + Type startAverage_; + + //- Current end index in sampleTimes + label endSampleTime_; + + //- Interpolated values from endSampleTime + Field<Type> endSampledValues_; + + //- If setAverage: end average value + Type endAverage_; + + +public: + + //- Runtime type information + TypeName("timeVaryingMappedFixedValue"); + + + // Constructors + + //- Construct from patch and internal field + timeVaryingMappedFixedValuePointPatchField + ( + const pointPatch&, + const DimensionedField<Type, pointMesh>& + ); + + //- Construct from patch, internal field and dictionary + timeVaryingMappedFixedValuePointPatchField + ( + const pointPatch&, + const DimensionedField<Type, pointMesh>&, + const dictionary& + ); + + //- Construct by mapping given patch field onto a new patch + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField<Type>&, + const pointPatch&, + const DimensionedField<Type, pointMesh>&, + const pointPatchFieldMapper& + ); + + //- Construct as copy + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField<Type>& + ); + + //- Construct and return a clone + virtual autoPtr<pointPatchField<Type> > clone() const + { + return autoPtr<pointPatchField<Type> > + ( + new timeVaryingMappedFixedValuePointPatchField<Type>(*this) + ); + } + + //- Construct as copy setting internal field reference + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField<Type>&, + const DimensionedField<Type, pointMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr<pointPatchField<Type> > clone + ( + const DimensionedField<Type, pointMesh>& iF + ) const + { + return autoPtr<pointPatchField<Type> > + ( + new timeVaryingMappedFixedValuePointPatchField<Type>(*this, iF) + ); + } + + + // Member functions + + // Utility functions + + //- Find boundary data inbetween current time and interpolate + void checkTable(); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "timeVaryingMappedFixedValuePointPatchField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C new file mode 100644 index 0000000000000000000000000000000000000000..ed847ab52d32588863e1b34e5c52c57f823018f1 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C @@ -0,0 +1,43 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "timeVaryingMappedFixedValuePointPatchFields.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makePointPatchFields(timeVaryingMappedFixedValue); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H new file mode 100644 index 0000000000000000000000000000000000000000..ad06b48066912a0505c9c6256ad3ff15c08925bb --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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/>. + +InClass + Foam::timeVaryingMappedFixedValuePointPatchFields + +Description + +SourceFiles + timeVaryingMappedFixedValuePointPatchFields.C + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingMappedFixedValuePointPatchFields_H +#define timeVaryingMappedFixedValuePointPatchFields_H + +#include "timeVaryingMappedFixedValuePointPatchField.H" +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchFieldTypedefs(timeVaryingMappedFixedValue); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 5912ea12e423618bc998e964cc4061d4be829be4..3329c2b28d1115a290a7cc0b4b11537fadc7cd11 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -158,6 +158,7 @@ triSurface/triangleFuncs/triangleFuncs.C triSurface/surfaceFeatures/surfaceFeatures.C triSurface/triSurfaceTools/triSurfaceTools.C triSurface/triSurfaceTools/geompack/geompack.C +triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C twoDPointCorrector/twoDPointCorrector.C diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C new file mode 100644 index 0000000000000000000000000000000000000000..c60be140affddc6c55d851763876fc3fad2da055 --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C @@ -0,0 +1,322 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "pointToPointPlanarInterpolation.H" +#include "boundBox.H" +#include "Random.H" +#include "vector2D.H" +#include "triSurface.H" +#include "triSurfaceTools.H" +#include "OFstream.H" +#include "Time.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::pointToPointPlanarInterpolation, 0); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::coordinateSystem +Foam::pointToPointPlanarInterpolation::calcCoordinateSystem +( + const pointField& points +) const +{ + if (points.size() < 3) + { + FatalErrorIn + ( + "pointToPointPlanarInterpolation::calcCoordinateSystem" + "(const pointField&)" + ) << "Only " << points.size() << " provided." << nl + << "Need at least three non-colinear points" + << " to be able to interpolate." + << exit(FatalError); + } + + const point& p0 = points[0]; + + // Find furthest away point + vector e1; + label index1 = -1; + scalar maxDist = -GREAT; + + for (label i = 1; i < points.size(); i++) + { + const vector d = points[i] - p0; + scalar magD = mag(d); + + if (magD > maxDist) + { + e1 = d/magD; + index1 = i; + maxDist = magD; + } + } + // Find point that is furthest away from line p0-p1 + const point& p1 = points[index1]; + + label index2 = -1; + maxDist = -GREAT; + for (label i = 1; i < points.size(); i++) + { + if (i != index1) + { + const point& p2 = points[i]; + vector e2(p2 - p0); + e2 -= (e2&e1)*e1; + scalar magE2 = mag(e2); + + if (magE2 > maxDist) + { + index2 = i; + maxDist = magE2; + } + } + } + if (index2 == -1) + { + FatalErrorIn + ( + "pointToPointPlanarInterpolation::calcCoordinateSystem" + "(const pointField&)" + ) << "Cannot find points that make valid normal." << nl + << "Have so far points " << p0 << " and " << p1 + << "Need at least three points which are not in a line." + << exit(FatalError); + } + + vector n = e1^(points[index2]-p0); + n /= mag(n); + + if (debug) + { + Info<< "pointToPointPlanarInterpolation::calcCoordinateSystem :" + << " Used points " << p0 << ' ' << points[index1] + << ' ' << points[index2] + << " to define coordinate system with normal " << n << endl; + } + + return coordinateSystem + ( + "reference", + p0, // origin + n, // normal + e1 // 0-axis + ); +} + + +void Foam::pointToPointPlanarInterpolation::calcWeights +( + const pointField& sourcePoints, + const pointField& destPoints +) +{ + tmp<vectorField> tlocalVertices + ( + referenceCS_.localPosition(sourcePoints) + ); + vectorField& localVertices = tlocalVertices(); + + const boundBox bb(localVertices, true); + const point bbMid(bb.midpoint()); + + if (debug) + { + Info<< "pointToPointPlanarInterpolation::readData :" + << " Perturbing points with " << perturb_ + << " fraction of a random position inside " << bb + << " to break any ties on regular meshes." + << nl << endl; + } + + Random rndGen(123456); + forAll(localVertices, i) + { + localVertices[i] += + perturb_ + *(rndGen.position(bb.min(), bb.max())-bbMid); + } + + // Determine triangulation + List<vector2D> localVertices2D(localVertices.size()); + forAll(localVertices, i) + { + localVertices2D[i][0] = localVertices[i][0]; + localVertices2D[i][1] = localVertices[i][1]; + } + + triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); + + tmp<pointField> tlocalFaceCentres + ( + referenceCS_.localPosition + ( + destPoints + ) + ); + const pointField& localFaceCentres = tlocalFaceCentres(); + + if (debug) + { + Pout<< "pointToPointPlanarInterpolation::readData :" + <<" Dumping triangulated surface to triangulation.stl" << endl; + s.write("triangulation.stl"); + + OFstream str("localFaceCentres.obj"); + Pout<< "readSamplePoints :" + << " Dumping face centres to " << str.name() << endl; + + forAll(localFaceCentres, i) + { + const point& p = localFaceCentres[i]; + str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; + } + } + + // Determine interpolation onto face centres. + triSurfaceTools::calcInterpolationWeights + ( + s, + localFaceCentres, // points to interpolate to + nearestVertex_, + nearestVertexWeight_ + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation +( + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb +) +: + perturb_(perturb), + referenceCS_(calcCoordinateSystem(sourcePoints)) + +{ + calcWeights(sourcePoints, destPoints); +} + + +Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation +( + const coordinateSystem& referenceCS, + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb +) +: + perturb_(perturb), + referenceCS_(referenceCS) +{ + calcWeights(sourcePoints, destPoints); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::wordList Foam::pointToPointPlanarInterpolation::timeNames +( + const instantList& times +) +{ + wordList names(times.size()); + + forAll(times, i) + { + names[i] = times[i].name(); + } + return names; +} + + +bool Foam::pointToPointPlanarInterpolation::findTime +( + const instantList& times, + const label startSampleTime, + const scalar timeVal, + label& lo, + label& hi +) +{ + lo = startSampleTime; + hi = -1; + + for (label i = startSampleTime+1; i < times.size(); i++) + { + if (times[i].value() > timeVal) + { + break; + } + else + { + lo = i; + } + } + + if (lo == -1) + { + //FatalErrorIn("findTime(..)") + // << "Cannot find starting sampling values for current time " + // << timeVal << nl + // << "Have sampling values for times " + // << timeNames(times) << nl + // << exit(FatalError); + return false; + } + + if (lo < times.size()-1) + { + hi = lo+1; + } + + + if (debug) + { + if (hi == -1) + { + Pout<< "findTime : Found time " << timeVal << " after" + << " index:" << lo << " time:" << times[lo].value() + << endl; + } + else + { + Pout<< "findTime : Found time " << timeVal << " inbetween" + << " index:" << lo << " time:" << times[lo].value() + << " and index:" << hi << " time:" << times[hi].value() + << endl; + } + } + return true; +} + + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H new file mode 100644 index 0000000000000000000000000000000000000000..daec2e46880a4aad4b79988e4bce3282f052f685 --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H @@ -0,0 +1,164 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::pointToPointPlanarInterpolation + +Description + Interpolates between two sets of unstructured points using 2D Delaunay + triangulation. Used in e.g. timeVaryingMapped bcs. + +SourceFiles + pointToPointPlanarInterpolation.C + +\*---------------------------------------------------------------------------*/ + +#ifndef pointToPointPlanarInterpolation_H +#define pointToPointPlanarInterpolation_H + +#include "FixedList.H" +#include "coordinateSystem.H" +#include "instantList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pointToPointPlanarInterpolation Declaration +\*---------------------------------------------------------------------------*/ + +class pointToPointPlanarInterpolation +{ + // Private data + + //- Perturbation factor + const scalar perturb_; + + //- Coordinate system + coordinateSystem referenceCS_; + + //- Current interpolation addressing to face centres of underlying + // patch + List<FixedList<label, 3> > nearestVertex_; + + //- Current interpolation factors to face centres of underlying + // patch + List<FixedList<scalar, 3> > nearestVertexWeight_; + + // Private Member Functions + + //- Calculate a local coordinate system from set of points + coordinateSystem calcCoordinateSystem(const pointField&) const; + + //- Calculate addressing and weights + void calcWeights + ( + const pointField& sourcePoints, + const pointField& destPoints + ); + +public: + + // Declare name of the class and its debug switch + ClassName("pointToPointPlanarInterpolation"); + + + // Constructors + + //- Construct from 3D locations. Determines local coordinate system + // from sourcePoints and maps onto that. + pointToPointPlanarInterpolation + ( + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb + ); + + //- Construct from coordinate system and locations. + pointToPointPlanarInterpolation + ( + const coordinateSystem& referenceCS, + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb + ); + + + // Member Functions + + //- Return the coordinateSystem + const coordinateSystem& referenceCS() const + { + return referenceCS_; + } + + // patch + const List<FixedList<label, 3> >& nearestVertex() const + { + return nearestVertex_; + } + + //- Current interpolation factors to face centres of underlying + // patch + const List<FixedList<scalar, 3> >& nearestVertexWeight() const + { + return nearestVertexWeight_; + } + + //- Helper: extract words of times + static wordList timeNames(const instantList&); + + //- Helper: find time. Return true if succesful. + static bool findTime + ( + const instantList& times, + const label startSampleTime, + const scalar timeVal, + label& lo, + label& hi + ); + + //- Interpolate from field on source points to dest points + template<class Type> + tmp<Field<Type> > interpolate(const Field<Type>& sourceFld) const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "pointToPointPlanarInterpolationTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..9159a03eddecdc0faef4261314ac203c1f536b6d --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C @@ -0,0 +1,71 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "pointToPointPlanarInterpolation.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +Foam::tmp<Foam::Field<Type> > Foam::pointToPointPlanarInterpolation::interpolate +( + const Field<Type>& sourceFld +) const +{ + tmp<Field<Type> > tfld(new Field<Type>(nearestVertex_.size())); + Field<Type>& fld = tfld(); + + forAll(fld, i) + { + const FixedList<label, 3>& verts = nearestVertex_[i]; + const FixedList<scalar, 3>& w = nearestVertexWeight_[i]; + + if (verts[2] == -1) + { + if (verts[1] == -1) + { + // Use vertex0 only + fld[i] = sourceFld[verts[0]]; + } + else + { + // Use vertex 0,1 + fld[i] = + w[0]*sourceFld[verts[0]] + + w[1]*sourceFld[verts[1]]; + } + } + else + { + fld[i] = + w[0]*sourceFld[verts[0]] + + w[1]*sourceFld[verts[1]] + + w[2]*sourceFld[verts[2]]; + } + } + return tfld; +} + + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H index 2328aa8284fc5478c5358ca631834c32c1832a24..ede548fd8349bf34835c7d6db3dec111345f33d2 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.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 @@ -25,7 +25,7 @@ Class Foam::triSurfaceTools Description - A collection of tools for triSurfaceMesh + A collection of tools for triSurface. SourceFiles triSurfaceTools.C