diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict index ba3dee7d..fbe33f18 100644 --- a/applications/utilities/postProcessing/sampling/sample/sampleDict +++ b/applications/utilities/postProcessing/sampling/sample/sampleDict @@ -141,6 +141,7 @@ sets patchSeed { type patchSeed; + axis xyz; patches (".*Wall.*"); // Number of points to seed. Divided amongst all processors according // to fraction of patches they hold. diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C index 6a6e6e0e..5a293d79 100644 --- a/applications/utilities/preProcessing/mapFields/mapFields.C +++ b/applications/utilities/preProcessing/mapFields/mapFields.C @@ -41,7 +41,8 @@ void mapConsistentMesh ( const fvMesh& meshSource, const fvMesh& meshTarget, - const meshToMesh::interpolationMethod& mapMethod, + const word& mapMethod, + const word& AMIMapMethod, const bool subtract, const HashSet& selectedFields, const bool noLagrangian @@ -50,7 +51,7 @@ void mapConsistentMesh Info<< nl << "Consistently creating and mapping fields for time " << meshSource.time().timeName() << nl << endl; - meshToMesh interp(meshSource, meshTarget, mapMethod); + meshToMesh interp(meshSource, meshTarget, mapMethod, AMIMapMethod); if (subtract) { @@ -79,7 +80,8 @@ void mapSubMesh const fvMesh& meshTarget, const HashTable& patchMap, const wordList& cuttingPatches, - const meshToMesh::interpolationMethod& mapMethod, + const word& mapMethod, + const word& AMIMapMethod, const bool subtract, const HashSet& selectedFields, const bool noLagrangian @@ -93,6 +95,7 @@ void mapSubMesh meshSource, meshTarget, mapMethod, + AMIMapMethod, patchMap, cuttingPatches ); @@ -186,6 +189,12 @@ int main(int argc, char *argv[]) "word", "specify the mapping method (direct|mapNearest|cellVolumeWeight)" ); + argList::addOption + ( + "patchMapMethod", + "word", + "specify the patch mapping method (direct|mapNearest|faceAreaWeight)" + ); argList::addBoolOption ( "subtract", @@ -231,17 +240,50 @@ int main(int argc, char *argv[]) const bool consistent = args.optionFound("consistent"); - meshToMesh::interpolationMethod mapMethod = - meshToMesh::imCellVolumeWeight; - if (args.optionFound("mapMethod")) + word mapMethod = meshToMesh::interpolationMethodNames_ + [ + meshToMesh::imCellVolumeWeight + ]; + + if (args.optionReadIfPresent("mapMethod", mapMethod)) { - mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]]; - - Info<< "Mapping method: " - << meshToMesh::interpolationMethodNames_[mapMethod] << endl; + Info<< "Mapping method: " << mapMethod << endl; } + + word patchMapMethod; + if (meshToMesh::interpolationMethodNames_.found(mapMethod)) + { + // Lookup corresponding AMI method + meshToMesh::interpolationMethod method = + meshToMesh::interpolationMethodNames_[mapMethod]; + + patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodToWord + ( + meshToMesh::interpolationMethodAMI(method) + ); + } + + // Optionally override + if (args.optionFound("patchMapMethod")) + { + patchMapMethod = args["patchMapMethod"]; + + Info<< "Patch mapping method: " << patchMapMethod << endl; + } + + + if (patchMapMethod.empty()) + { + FatalErrorIn(args.executable()) + << "No valid patchMapMethod for method " << mapMethod + << ". Please supply one through the 'patchMapMethod' option" + << exit(FatalError); + } + + + const bool subtract = args.optionFound("subtract"); if (subtract) { @@ -314,6 +356,7 @@ int main(int argc, char *argv[]) meshSource, meshTarget, mapMethod, + patchMapMethod, subtract, selectedFields, noLagrangian @@ -328,6 +371,7 @@ int main(int argc, char *argv[]) patchMap, addProcessorPatches(meshTarget, cuttingPatches), mapMethod, + patchMapMethod, subtract, selectedFields, noLagrangian diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 668a5953..945e727b 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1179,7 +1179,9 @@ Foam::tmp Foam::polyMesh::movePoints if (debug && moveError) { - write(); + // Write mesh to ease debugging. Note we want to avoid calling + // e.g. fvMesh::write since meshPhi not yet complete. + polyMesh::write(); } return sweptVols; diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C index f93817c3..48c3af19 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-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,6 +72,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), + mapMethod_(ptf.mapMethod_), mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), @@ -102,6 +103,14 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(iF.name()), setAverage_(readBool(dict.lookup("setAverage"))), perturb_(dict.lookupOrDefault("perturb", 1e-5)), + mapMethod_ + ( + dict.lookupOrDefault + ( + "mapMethod", + "planarInterpolation" + ) + ), mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), @@ -112,6 +121,27 @@ timeVaryingMappedFixedValueFvPatchField endAverage_(pTraits::zero), offset_(DataEntry::New("offset", dict)) { + if + ( + mapMethod_ != "planarInterpolation" + && mapMethod_ != "nearest" + ) + { + FatalIOErrorIn + ( + "timeVaryingMappedFixedValueFvPatchField::\n" + "timeVaryingMappedFixedValueFvPatchField\n" + "(\n" + " const fvPatch&\n" + " const DimensionedField&\n" + " const dictionary&\n" + ")\n", + dict + ) << "mapMethod should be one of 'planarInterpolation'" + << ", 'nearest'" << exit(FatalIOError); + } + + dict.readIfPresent("fieldTableName", fieldTableName_); if (dict.found("value")) @@ -140,6 +170,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), + mapMethod_(ptf.mapMethod_), mapperPtr_(NULL), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), @@ -169,6 +200,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), + mapMethod_(ptf.mapMethod_), mapperPtr_(NULL), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), @@ -258,6 +290,14 @@ void timeVaryingMappedFixedValueFvPatchField::checkTable() << samplePointsFile << endl; } + + // tbd: run-time selection + bool nearestOnly = + ( + !mapMethod_.empty() + && mapMethod_ != "planarInterpolation" + ); + // Allocate the interpolator mapperPtr_.reset ( @@ -265,7 +305,8 @@ void timeVaryingMappedFixedValueFvPatchField::checkTable() ( samplePoints, this->patch().patch().faceCentres(), - perturb_ + perturb_, + nearestOnly ) ); @@ -560,6 +601,18 @@ void timeVaryingMappedFixedValueFvPatchField::write(Ostream& os) const << token::END_STATEMENT << nl; } + if + ( + ( + !mapMethod_.empty() + && mapMethod_ != "planarInterpolation" + ) + ) + { + os.writeKeyword("mapMethod") << mapMethod_ + << token::END_STATEMENT << nl; + } + offset_->writeData(os); this->writeEntry("value", os); diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H index aec6c322..8237cdde 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-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -33,12 +33,13 @@ Description constant/boundaryData/\ where: - points : pointField with locations - ddd : supplied values at time ddd - The points should be more or less on a plane since they get triangulated - in 2-D. + The default mode of operation (mapMethod planarInterpolation) is + to project the points onto a plane (constructed from the first threee + points) and construct a 2D triangulation and finds for the face centres + the triangle it is in and the weights to the 3 vertices. - At startup, this condition generates the triangulation and performs a - linear interpolation (triangle it is in and weights to the 3 vertices) - for every face centre. + The optional mapMethod nearest will avoid all projection and + triangulation and just use the value at the nearest vertex. Values are interpolated linearly between times. @@ -49,6 +50,7 @@ Description setAverage | flag to activate setting of average value | yes | perturb | perturb points for regular geometries | no | 1e-5 fieldTableName | alternative field name to sample | no| this field name + mapMethod | type of mapping | no | planarInterpolation \endtable /verbatim @@ -61,10 +63,6 @@ Description } /endverbatim -Note - Switch on debug flag to have it dump the triangulation (in transformed - space) and transform face centres. - SeeAlso Foam::fixedValueFvPatchField @@ -107,7 +105,10 @@ class timeVaryingMappedFixedValueFvPatchField //- Fraction of perturbation (fraction of bounding box) to add scalar perturb_; - //- 2D interpolation + //- Interpolation scheme to use + word mapMethod_; + + //- 2D interpolation (for 'planarInterpolation' mapMethod) autoPtr mapperPtr_; //- List of boundaryData time directories diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C index bbbc75f3..54725088 100644 --- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -42,6 +42,7 @@ timeVaryingMappedFixedValuePointPatchField fieldTableName_(iF.name()), setAverage_(false), perturb_(0), + mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -68,7 +69,8 @@ timeVaryingMappedFixedValuePointPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - mapperPtr_(ptf.mapperPtr_), + mapMethod_(ptf.mapMethod_), + mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -99,6 +101,14 @@ timeVaryingMappedFixedValuePointPatchField fieldTableName_(iF.name()), setAverage_(readBool(dict.lookup("setAverage"))), perturb_(dict.lookupOrDefault("perturb", 1e-5)), + mapMethod_ + ( + dict.lookupOrDefault + ( + "mapMethod", + "planarInterpolation" + ) + ), mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), @@ -146,6 +156,7 @@ timeVaryingMappedFixedValuePointPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), + mapMethod_(ptf.mapMethod_), mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), @@ -176,6 +187,7 @@ timeVaryingMappedFixedValuePointPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), + mapMethod_(ptf.mapMethod_), mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), @@ -290,13 +302,22 @@ void Foam::timeVaryingMappedFixedValuePointPatchField::checkTable() ) ); + // tbd: run-time selection + bool nearestOnly = + ( + !mapMethod_.empty() + && mapMethod_ != "planarInterpolation" + ); + + // Allocate the interpolator mapperPtr_.reset ( new pointToPointPlanarInterpolation ( samplePoints, meshPts, - perturb_ + perturb_, + nearestOnly ) ); @@ -593,6 +614,18 @@ void Foam::timeVaryingMappedFixedValuePointPatchField::write << token::END_STATEMENT << nl; } + if + ( + ( + !mapMethod_.empty() + && mapMethod_ != "planarInterpolation" + ) + ) + { + os.writeKeyword("mapMethod") << mapMethod_ + << token::END_STATEMENT << nl; + } + if (offset_.valid()) { offset_->writeData(os); diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H index ed87fa4c..012ac08d 100644 --- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,7 +68,10 @@ class timeVaryingMappedFixedValuePointPatchField //- Fraction of perturbation (fraction of bounding box) to add scalar perturb_; - //- 2D interpolation + //- Interpolation scheme to use + word mapMethod_; + + //- 2D interpolation (for 'planarInterpolation' mapMethod) autoPtr mapperPtr_; //- List of boundaryData time directories diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index 08c8f7e5..4999a16d 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -2066,6 +2066,15 @@ void Foam::autoLayerDriver::setupLayerInfoTruncation } } nPatchPointLayers = patchNLayers; + + // Set any unset patch face layers + forAll(nPatchFaceLayers, patchFaceI) + { + if (nPatchFaceLayers[patchFaceI] == -1) + { + nPatchFaceLayers[patchFaceI] = 0; + } + } } else { diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index ea79e069..58fb989e 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -536,66 +536,13 @@ void Foam::AMIInterpolation::agglomerate } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - template -Foam::AMIInterpolation::AMIInterpolation +void Foam::AMIInterpolation::constructFromSurface ( const SourcePatch& srcPatch, const TargetPatch& tgtPatch, - const faceAreaIntersect::triangulationMode& triMode, - const bool requireMatch, - const interpolationMethod& method, - const scalar lowWeightCorrection, - const bool reverseTarget + const autoPtr& surfPtr ) -: - method_(method), - reverseTarget_(reverseTarget), - requireMatch_(requireMatch), - singlePatchProc_(-999), - lowWeightCorrection_(lowWeightCorrection), - srcAddress_(), - srcWeights_(), - srcWeightsSum_(), - tgtAddress_(), - tgtWeights_(), - tgtWeightsSum_(), - triMode_(triMode), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) -{ - update(srcPatch, tgtPatch); -} - - -template -Foam::AMIInterpolation::AMIInterpolation -( - const SourcePatch& srcPatch, - const TargetPatch& tgtPatch, - const autoPtr& surfPtr, - const faceAreaIntersect::triangulationMode& triMode, - const bool requireMatch, - const interpolationMethod& method, - const scalar lowWeightCorrection, - const bool reverseTarget -) -: - method_(method), - reverseTarget_(reverseTarget), - requireMatch_(requireMatch), - singlePatchProc_(-999), - lowWeightCorrection_(lowWeightCorrection), - srcAddress_(), - srcWeights_(), - srcWeightsSum_(), - tgtAddress_(), - tgtWeights_(), - tgtWeightsSum_(), - triMode_(triMode), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) { if (surfPtr.valid()) { @@ -658,6 +605,134 @@ Foam::AMIInterpolation::AMIInterpolation } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const interpolationMethod& method, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(interpolationMethodToWord(method)), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + update(srcPatch, tgtPatch); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const word& methodName, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(methodName), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + update(srcPatch, tgtPatch); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const interpolationMethod& method, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(interpolationMethodToWord(method)), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + constructFromSurface(srcPatch, tgtPatch, surfPtr); +} + + +template +Foam::AMIInterpolation::AMIInterpolation +( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch, + const word& methodName, + const scalar lowWeightCorrection, + const bool reverseTarget +) +: + methodName_(methodName), + reverseTarget_(reverseTarget), + requireMatch_(requireMatch), + singlePatchProc_(-999), + lowWeightCorrection_(lowWeightCorrection), + srcAddress_(), + srcWeights_(), + srcWeightsSum_(), + tgtAddress_(), + tgtWeights_(), + tgtWeightsSum_(), + triMode_(triMode), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + constructFromSurface(srcPatch, tgtPatch, surfPtr); +} + + template Foam::AMIInterpolation::AMIInterpolation ( @@ -666,7 +741,7 @@ Foam::AMIInterpolation::AMIInterpolation const labelList& targetRestrictAddressing ) : - method_(fineAMI.method_), + methodName_(fineAMI.methodName_), reverseTarget_(fineAMI.reverseTarget_), requireMatch_(fineAMI.requireMatch_), singlePatchProc_(fineAMI.singlePatchProc_), @@ -883,7 +958,7 @@ void Foam::AMIInterpolation::update ( AMIMethod::New ( - interpolationMethodToWord(method_), + methodName_, srcPatch, newTgtPatch, srcMagSf_, @@ -1000,7 +1075,7 @@ void Foam::AMIInterpolation::update ( AMIMethod::New ( - interpolationMethodToWord(method_), + methodName_, srcPatch, tgtPatch, srcMagSf_, diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index 7dbacd9e..27d077f1 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -110,7 +110,7 @@ private: // Private data //- Interpolation method - interpolationMethod method_; + const word methodName_; //- Flag to indicate that the two patches are co-directional and // that the orientation of the target patch should be reversed @@ -250,7 +250,7 @@ private: ); - // Constructor helper + // Constructor helpers static void agglomerate ( @@ -269,6 +269,12 @@ private: autoPtr& tgtMap ); + void constructFromSurface + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surfPtr + ); public: @@ -286,6 +292,19 @@ public: const bool reverseTarget = false ); + //- Construct from components + AMIInterpolation + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch = true, + const word& methodName = + interpolationMethodToWord(imFaceAreaWeight), + const scalar lowWeightCorrection = -1, + const bool reverseTarget = false + ); + //- Construct from components, with projection surface AMIInterpolation ( @@ -299,6 +318,20 @@ public: const bool reverseTarget = false ); + //- Construct from components, with projection surface + AMIInterpolation + ( + const SourcePatch& srcPatch, + const TargetPatch& tgtPatch, + const autoPtr& surf, + const faceAreaIntersect::triangulationMode& triMode, + const bool requireMatch = true, + const word& methodName = + interpolationMethodToWord(imFaceAreaWeight), + const scalar lowWeightCorrection = -1, + const bool reverseTarget = false + ); + //- Construct from agglomeration of AMIInterpolation. Agglomeration // passed in as new coarse size and addressing from fine from coarse AMIInterpolation diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options index 1b22dab2..1713152e 100644 --- a/src/meshTools/Make/options +++ b/src/meshTools/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/triSurface/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ -I$(LIB_SRC)/fileFormats/lnInclude LIB_LIBS = \ -ltriSurface \ + -lsurfMesh \ -lfileFormats diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C index 1f7de7c5..5d69ba1f 100644 --- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,8 +29,9 @@ License #include "vector2D.H" #include "triSurface.H" #include "triSurfaceTools.H" -#include "OFstream.H" +#include "OBJstream.H" #include "Time.H" +#include "matchPoints.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -139,114 +140,172 @@ void Foam::pointToPointPlanarInterpolation::calcWeights const pointField& destPoints ) { - tmp tlocalVertices - ( - referenceCS_.localPosition(sourcePoints) - ); - vectorField& localVertices = tlocalVertices(); - - const boundBox bb(localVertices, true); - const point bbMid(bb.midpoint()); - - if (debug) + if (nearestOnly_) { - 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 localVertices2D(localVertices.size()); - forAll(localVertices, i) - { - localVertices2D[i][0] = localVertices[i][0]; - localVertices2D[i][1] = localVertices[i][1]; - } - - triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); - - tmp tlocalFaceCentres - ( - referenceCS_.localPosition + labelList destToSource; + bool fullMatch = matchPoints ( - destPoints - ) - ); - const pointField& localFaceCentres = tlocalFaceCentres(); + destPoints, + sourcePoints, + scalarField(destPoints.size(), GREAT), + true, // verbose + destToSource + ); - 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) + if (!fullMatch) { - const point& p = localFaceCentres[i]; - str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; + FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)") + << "Did not find a corresponding sourcePoint for every face" + << " centre" << exit(FatalError); + } + + nearestVertex_.setSize(destPoints.size()); + nearestVertexWeight_.setSize(destPoints.size()); + forAll(nearestVertex_, i) + { + nearestVertex_[i][0] = destToSource[i]; + nearestVertex_[i][1] = -1; + nearestVertex_[i][2] = -1; + + nearestVertexWeight_[i][0] = 1.0; + nearestVertexWeight_[i][1] = 0.0; + nearestVertexWeight_[i][2] = 0.0; + } + + if (debug) + { + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + + Pout<< "For location " << destPoints[i] + << " sampling vertex " << v0 + << " at:" << sourcePoints[v0] + << " distance:" << mag(sourcePoints[v0]-destPoints[i]) + << endl; + } + + OBJstream str("destToSource.obj"); + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + << " Dumping lines from face centres to original points to " + << str.name() << endl; + + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + str.write(linePointRef(destPoints[i], sourcePoints[v0])); + } } } - - // Determine interpolation onto face centres. - triSurfaceTools::calcInterpolationWeights - ( - s, - localFaceCentres, // points to interpolate to - nearestVertex_, - nearestVertexWeight_ - ); - - if (debug) + else { - forAll(sourcePoints, i) + tmp tlocalVertices + ( + referenceCS_.localPosition(sourcePoints) + ); + vectorField& localVertices = tlocalVertices(); + + const boundBox bb(localVertices, true); + const point bbMid(bb.midpoint()); + + if (debug) { - Pout<< "source:" << i << " at:" << sourcePoints[i] - << " 2d:" << localVertices[i] - << endl; + Info<< "pointToPointPlanarInterpolation::calcWeights :" + << " Perturbing points with " << perturb_ + << " fraction of a random position inside " << bb + << " to break any ties on regular meshes." + << nl << endl; } - - forAll(destPoints, i) + Random rndGen(123456); + forAll(localVertices, i) { - label v0 = nearestVertex_[i][0]; - label v1 = nearestVertex_[i][1]; - label v2 = nearestVertex_[i][2]; + localVertices[i] += + perturb_ + *(rndGen.position(bb.min(), bb.max())-bbMid); + } - Pout<< "For location " << destPoints[i] - << " 2d:" << localFaceCentres[i] - << " sampling vertices" << nl - << " " << v0 - << " at:" << sourcePoints[v0] - << " weight:" << nearestVertexWeight_[i][0] << nl; + // Determine triangulation + List localVertices2D(localVertices.size()); + forAll(localVertices, i) + { + localVertices2D[i][0] = localVertices[i][0]; + localVertices2D[i][1] = localVertices[i][1]; + } - if (v1 != -1) + triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); + + tmp tlocalFaceCentres + ( + referenceCS_.localPosition + ( + destPoints + ) + ); + const pointField& localFaceCentres = tlocalFaceCentres(); + + if (debug) + { + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + <<" Dumping triangulated surface to triangulation.stl" << endl; + s.write("triangulation.stl"); + + OBJstream str("localFaceCentres.obj"); + Pout<< "pointToPointPlanarInterpolation::calcWeights :" + << " Dumping face centres to " << str.name() << endl; + + forAll(localFaceCentres, i) { - Pout<< " " << v1 - << " at:" << sourcePoints[v1] - << " weight:" << nearestVertexWeight_[i][1] << nl; + str.write(localFaceCentres[i]); } - if (v2 != -1) + } + + // Determine interpolation onto face centres. + triSurfaceTools::calcInterpolationWeights + ( + s, + localFaceCentres, // points to interpolate to + nearestVertex_, + nearestVertexWeight_ + ); + + if (debug) + { + forAll(sourcePoints, i) { - Pout<< " " << v2 - << " at:" << sourcePoints[v2] - << " weight:" << nearestVertexWeight_[i][2] << nl; + Pout<< "source:" << i << " at:" << sourcePoints[i] + << " 2d:" << localVertices[i] + << endl; } - Pout<< endl; + forAll(destPoints, i) + { + label v0 = nearestVertex_[i][0]; + label v1 = nearestVertex_[i][1]; + label v2 = nearestVertex_[i][2]; + + Pout<< "For location " << destPoints[i] + << " 2d:" << localFaceCentres[i] + << " sampling vertices" << nl + << " " << v0 + << " at:" << sourcePoints[v0] + << " weight:" << nearestVertexWeight_[i][0] << nl; + + if (v1 != -1) + { + Pout<< " " << v1 + << " at:" << sourcePoints[v1] + << " weight:" << nearestVertexWeight_[i][1] << nl; + } + if (v2 != -1) + { + Pout<< " " << v2 + << " at:" << sourcePoints[v2] + << " weight:" << nearestVertexWeight_[i][2] << nl; + } + + Pout<< endl; + } } } } @@ -258,13 +317,14 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation ( const pointField& sourcePoints, const pointField& destPoints, - const scalar perturb + const scalar perturb, + const bool nearestOnly ) : perturb_(perturb), + nearestOnly_(nearestOnly), referenceCS_(calcCoordinateSystem(sourcePoints)), nPoints_(sourcePoints.size()) - { calcWeights(sourcePoints, destPoints); } @@ -279,6 +339,7 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation ) : perturb_(perturb), + nearestOnly_(false), referenceCS_(referenceCS), nPoints_(sourcePoints.size()) { diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H index e3428d80..4683dcd4 100644 --- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -56,6 +56,9 @@ class pointToPointPlanarInterpolation //- Perturbation factor const scalar perturb_; + //- Whether to use nearest point only (avoids triangulation, projection) + const bool nearestOnly_; + //- Coordinate system coordinateSystem referenceCS_; @@ -91,12 +94,15 @@ public: // Constructors //- Construct from 3D locations. Determines local coordinate system - // from sourcePoints and maps onto that. + // from sourcePoints and maps onto that. If nearestOnly skips any + // local coordinate system and triangulation and uses nearest vertex + // only pointToPointPlanarInterpolation ( const pointField& sourcePoints, const pointField& destPoints, - const scalar perturb + const scalar perturb, + const bool nearestOnly = false ); //- Construct from coordinate system and locations. diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C index 2b2b7e97..91c998ce 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C +++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C @@ -119,6 +119,7 @@ void Foam::meshToMesh::normaliseWeights void Foam::meshToMesh::calcAddressing ( + const word& methodName, const polyMesh& src, const polyMesh& tgt ) @@ -127,7 +128,7 @@ void Foam::meshToMesh::calcAddressing ( meshToMeshMethod::New ( - interpolationMethodNames_[method_], + methodName, src, tgt ) @@ -150,11 +151,11 @@ void Foam::meshToMesh::calcAddressing } -void Foam::meshToMesh::calculate() +void Foam::meshToMesh::calculate(const word& methodName) { Info<< "Creating mesh-to-mesh addressing for " << srcRegion_.name() << " and " << tgtRegion_.name() << " regions using " - << interpolationMethodNames_[method_] << endl; + << methodName << endl; singleMeshProc_ = calcDistribution(srcRegion_, tgtRegion_); @@ -241,7 +242,7 @@ void Foam::meshToMesh::calculate() } } - calcAddressing(srcRegion_, newTgt); + calcAddressing(methodName, srcRegion_, newTgt); // per source cell the target cell address in newTgt mesh forAll(srcToTgtCellAddr_, i) @@ -320,7 +321,7 @@ void Foam::meshToMesh::calculate() } else { - calcAddressing(srcRegion_, tgtRegion_); + calcAddressing(methodName, srcRegion_, tgtRegion_); normaliseWeights ( @@ -342,12 +343,9 @@ void Foam::meshToMesh::calculate() Foam::AMIPatchToPatchInterpolation::interpolationMethod -Foam::meshToMesh::interpolationMethodAMI -( - const interpolationMethod method -) const +Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method) { - switch (method_) + switch (method) { case imDirect: { @@ -374,7 +372,7 @@ Foam::meshToMesh::interpolationMethodAMI "const interpolationMethod method" ") const" ) - << "Unhandled enumeration " << method_ + << "Unhandled enumeration " << method << abort(FatalError); } } @@ -383,90 +381,66 @@ Foam::meshToMesh::interpolationMethodAMI } -const Foam::PtrList& -Foam::meshToMesh::patchAMIs() const +void Foam::meshToMesh::calculatePatchAMIs(const word& AMIMethodName) { - if (patchAMIs_.empty()) + if (!patchAMIs_.empty()) { - const word amiMethod = - AMIPatchToPatchInterpolation::interpolationMethodToWord - ( - interpolationMethodAMI(method_) - ); - - patchAMIs_.setSize(srcPatchID_.size()); - - forAll(srcPatchID_, i) - { - label srcPatchI = srcPatchID_[i]; - label tgtPatchI = tgtPatchID_[i]; - - const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI]; - const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI]; - - Info<< "Creating AMI between source patch " << srcPP.name() - << " and target patch " << tgtPP.name() - << " using " << amiMethod - << endl; - - Info<< incrIndent; - - patchAMIs_.set - ( - i, - new AMIPatchToPatchInterpolation - ( - srcPP, - tgtPP, - faceAreaIntersect::tmMesh, - false, - interpolationMethodAMI(method_), - -1, - true // flip target patch since patch normals are aligned - ) - ); - - Info<< decrIndent; - } + FatalErrorIn("meshToMesh::calculatePatchAMIs()") + << "patch AMI already calculated" + << exit(FatalError); } - return patchAMIs_; + patchAMIs_.setSize(srcPatchID_.size()); + + forAll(srcPatchID_, i) + { + label srcPatchI = srcPatchID_[i]; + label tgtPatchI = tgtPatchID_[i]; + + const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI]; + const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI]; + + Info<< "Creating AMI between source patch " << srcPP.name() + << " and target patch " << tgtPP.name() + << " using " << AMIMethodName + << endl; + + Info<< incrIndent; + + patchAMIs_.set + ( + i, + new AMIPatchToPatchInterpolation + ( + srcPP, + tgtPP, + faceAreaIntersect::tmMesh, + false, + AMIMethodName, + -1, + true // flip target patch since patch normals are aligned + ) + ); + + Info<< decrIndent; + } } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::meshToMesh::meshToMesh +void Foam::meshToMesh::constructNoCuttingPatches ( - const polyMesh& src, - const polyMesh& tgt, - const interpolationMethod& method, - bool interpAllPatches + const word& methodName, + const word& AMIMethodName, + const bool interpAllPatches ) -: - srcRegion_(src), - tgtRegion_(tgt), - srcPatchID_(), - tgtPatchID_(), - patchAMIs_(), - cuttingPatches_(), - srcToTgtCellAddr_(), - tgtToSrcCellAddr_(), - srcToTgtCellWght_(), - tgtToSrcCellWght_(), - method_(method), - V_(0.0), - singleMeshProc_(-1), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) { if (interpAllPatches) { - const polyBoundaryMesh& srcBM = src.boundaryMesh(); - const polyBoundaryMesh& tgtBM = tgt.boundaryMesh(); + const polyBoundaryMesh& srcBM = srcRegion_.boundaryMesh(); + const polyBoundaryMesh& tgtBM = tgtRegion_.boundaryMesh(); - DynamicList