diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index b58df5402e39ac354eadf4f075ca861156d595b9..67df5b61e054ef46b32c8daefe3e367023496248 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -22,7 +22,7 @@ License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Application - buoyantSimpleRadiationFoam + buoyantSimpleFoam Description Steady-state solver for buoyant, turbulent flow of compressible fluids, diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options index 06126f0bf86d8f3b6fe8545a230cb91b899dcff3..94dd40ce6b944099f73a4f5b643f4426d90fad1e 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options @@ -11,9 +11,7 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \ - -I$(LIB_SRC)/turbulenceModels \ -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \ -I$(LIB_SRC)/fvOptions/lnInclude \ -I$(LIB_SRC)/regionModels/regionModel/lnInclude @@ -25,6 +23,7 @@ EXE_LIBS = \ -lspecie \ -lcompressibleTurbulenceModel \ -lcompressibleRASModels \ + -lcompressibleLESModels \ -lmeshTools \ -lfiniteVolume \ -lradiationModels \ diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C index 7e45266d241298f134d0232c163651160fafadda..137339e631ee2d9c3153883b8f52430583e2fc70 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C @@ -150,13 +150,13 @@ int Foam::processorCyclicPolyPatch::tag() const referPatch() ); - if (cycPatch.owner()) + if (owner()) { - tag_ = Hash<word>()(cycPatch.name()); + tag_ = Hash<word>()(cycPatch.name()) % 32768u; } else { - tag_ = Hash<word>()(cycPatch.neighbPatch().name()); + tag_ = Hash<word>()(cycPatch.neighbPatch().name()) % 32768u; } if (tag_ == Pstream::msgType() || tag_ == -1) diff --git a/src/OpenFOAM/primitives/ops/ops.H b/src/OpenFOAM/primitives/ops/ops.H index cab91fcfb5eee81d71a4cb8280b8d21e0c25e1a0..427679118f2e2a0f179285f8855da314fac8cb6c 100644 --- a/src/OpenFOAM/primitives/ops/ops.H +++ b/src/OpenFOAM/primitives/ops/ops.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -175,7 +175,7 @@ Op(lessEq, x <= y) Op(greater, x > y) Op(greaterEq, x >= y) -weightedOp(multiply, weight * y) +weightedOp(multiply, (weight*y)) #undef Op diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C index 8bca126f9b503a10f0ff1488117626d4107f74a3..edcd775d3a9fb8e1f928a9bb108ea75606b94fb4 100644 --- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C @@ -156,6 +156,47 @@ timeVaryingMappedFixedValuePointPatchField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template<class Type> +void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::autoMap +( + const pointPatchFieldMapper& m +) +{ + fixedValuePointPatchField<Type>::autoMap(m); + if (startSampledValues_.size()) + { + startSampledValues_.autoMap(m); + endSampledValues_.autoMap(m); + } + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; +} + + +template<class Type> +void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::rmap +( + const pointPatchField<Type>& ptf, + const labelList& addr +) +{ + fixedValuePointPatchField<Type>::rmap(ptf, addr); + + const timeVaryingMappedFixedValuePointPatchField<Type>& tiptf = + refCast<const timeVaryingMappedFixedValuePointPatchField<Type> >(ptf); + + startSampledValues_.rmap(tiptf.startSampledValues_, addr); + endSampledValues_.rmap(tiptf.endSampledValues_, addr); + + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; +} + + template<class Type> void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable() { diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H index 108ec6786ee9b4a804ec0e5c1ebe0068f56ae538..942dc3e4a643960058160516e32533b59cc49937 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 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -167,6 +167,23 @@ public: void checkTable(); + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const pointPatchFieldMapper& + ); + + //- Reverse map the given PointPatchField onto + // this PointPatchField + virtual void rmap + ( + const pointPatchField<Type>&, + const labelList& + ); + + // Evaluation functions //- Update the coefficients associated with the patch field diff --git a/src/fvOptions/fvOptions/fvOption.C b/src/fvOptions/fvOptions/fvOption.C index 8cb72b8e996b90a89826b037d0781aa5f6b4da5b..1a46a99794a25dd36132dca04e79be3130856857 100644 --- a/src/fvOptions/fvOptions/fvOption.C +++ b/src/fvOptions/fvOptions/fvOption.C @@ -193,7 +193,8 @@ void Foam::fv::option::setCellSet() meshToMeshNew::interpolationMethodNames_.read ( dict_.lookup("interpolationMethod") - ) + ), + false // not interpolating patches ) ); } diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C index f535f93be40d8bc6ab9a944490f12afe22ec8ba8..d1b085f1574d0db0739ad138d8260641c57beedf 100644 --- a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C +++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C @@ -78,6 +78,9 @@ void Foam::fv::interRegionHeatTransferModel::setNbrModel() } firstIter_ = false; + + // set nbr model's nbr model to avoid construction order problems + nbrModel_->setNbrModel(); } diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C index e72f13d72663e3683a878f7267c116d2428be75e..9c8f08630f189deeda3eacee7cc520896c58f7c0 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C @@ -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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,23 +41,31 @@ void Foam::refinementFeatures::read fileName featFileName(dict.lookup("file")); - set - ( - featI, - new featureEdgeMesh + { + IOobject featObj ( - IOobject + featFileName, // name + io.time().constant(), // instance + "triSurface", // local + io.time(), // registry + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ); + + autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(featObj.filePath()); + + set + ( + featI, + new featureEdgeMesh ( - featFileName, // name - io.time().constant(), // instance - "triSurface", // local - io.time(), // registry - IOobject::MUST_READ, - IOobject::NO_WRITE, - false + featObj, + eMeshPtr->points(), + eMeshPtr->edges() ) - ) - ); + ); + } const featureEdgeMesh& eMesh = operator[](featI); diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index fc60cd2fbba5c8d30cac94380afff83e05c7a7c9..f62ea404d9ae2090ea452e6f3fbadcf4d4001ba1 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -857,7 +857,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights if (nFace) { - Info<< "AMI: Patch " << patchName << " weights min/max/average = " + IInfo<< "AMI: Patch " << patchName << " weights min/max/average = " << gMin(wghtSum) << ", " << gMax(wghtSum) << ", " << gAverage(wghtSum) << endl; @@ -1156,7 +1156,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation label srcSize = returnReduce(srcPatch.size(), sumOp<label>()); label tgtSize = returnReduce(tgtPatch.size(), sumOp<label>()); - Info<< "AMI: Creating addressing and weights between " + IInfo<< "AMI: Creating addressing and weights between " << srcSize << " source faces and " << tgtSize << " target faces" << endl; @@ -1191,7 +1191,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation label srcSize = returnReduce(srcPatch.size(), sumOp<label>()); label tgtSize = returnReduce(tgtPatch.size(), sumOp<label>()); - Info<< "AMI: Creating addressing and weights between " + IInfo<< "AMI: Creating addressing and weights between " << srcSize << " source faces and " << tgtSize << " target faces" << endl; @@ -1312,9 +1312,9 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation ( "AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation" "(" - " const AMIInterpolation<SourcePatch, TargetPatch>&, " - " const labelList&, " - " const labelList&" + "const AMIInterpolation<SourcePatch, TargetPatch>&, " + "const labelList&, " + "const labelList&" ")" ) << "Size mismatch." << nl << "Source patch size:" << fineAMI.srcAddress().size() << nl @@ -1579,7 +1579,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update if (debug) { Info<< "AMIInterpolation : Constructed addressing and weights" << nl - << " triMode :" << triMode_ << nl + << " triMode :" + << faceAreaIntersect::triangulationModeNames_[triMode_] << nl << " singlePatchProc:" << singlePatchProc_ << nl << " srcMagSf :" << gSum(srcMagSf_) << nl << " tgtMagSf :" << gSum(tgtMagSf_) << nl diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index 0b7633035536627ff5134a06e76589f5fd63eb89..9f94fa41169965fbd12bece06b9c0afc0edb6563 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -469,7 +469,7 @@ public: ) const; - //- Interpolate from target to source with supplied binary op + //- Interpolate from target to source with supplied op template<class Type, class CombineOp> tmp<Field<Type> > interpolateToSource ( @@ -477,8 +477,7 @@ public: const CombineOp& cop ) const; - //- Interpolate from target tmp field to source with supplied - // binary op + //- Interpolate from target tmp field to source with supplied op template<class Type, class CombineOp> tmp<Field<Type> > interpolateToSource ( diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C index 4135cfe04191b4ffc2e2135c41dbf5b802a6b5dc..45d1eb515b519035e4c9f733709636414890899c 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C @@ -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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -27,6 +27,19 @@ License // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // +namespace Foam +{ + template<> + const char* NamedEnum<faceAreaIntersect::triangulationMode, 2>::names[] = + { + "fan", + "mesh" + }; +} + +const Foam::NamedEnum<Foam::faceAreaIntersect::triangulationMode, 2> + Foam::faceAreaIntersect::triangulationModeNames_; + Foam::scalar Foam::faceAreaIntersect::tol = 1e-6; // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H index bfbe7c08a2e2a1d31cbec044e7a01d12eea68399..a8dd739151f9a2a1756b6a3978b2952ae07ec370 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,6 +41,7 @@ SourceFiles #include "FixedList.H" #include "plane.H" #include "face.H" +#include "NamedEnum.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -63,6 +64,8 @@ public: tmMesh }; + static const NamedEnum<triangulationMode, 2> triangulationModeNames_; + private: diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/force/force.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/force/force.H index 24e921da73d122af4ceefe69649c1f75513e7ef0..59389c387bf2cedd499780826436e4252b438f8d 100644 --- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/force/force.H +++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/force/force.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -108,7 +108,7 @@ public: ( const surfaceFilmModel& owner, const dictionary& dict, - const word& mdoelType + const word& modelType ); diff --git a/src/sampling/Make/files b/src/sampling/Make/files index 747097a781c5c26d52c5ae82c95b484bed4cd0ae..309cf09d7ab39e79b984b4c0e7e8cce10f873b09 100644 --- a/src/sampling/Make/files +++ b/src/sampling/Make/files @@ -60,6 +60,9 @@ $(meshToMesh)/calculateMeshToMeshAddressing.C $(meshToMesh)/calculateMeshToMeshWeights.C meshToMeshNew = meshToMeshInterpolation/meshToMeshNew +$(meshToMeshNew)/calcDirect.C +$(meshToMeshNew)/calcMapNearest.C +$(meshToMeshNew)/calcCellVolumeWeight.C $(meshToMeshNew)/meshToMeshNew.C $(meshToMeshNew)/meshToMeshNewParallelOps.C diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcCellVolumeWeight.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcCellVolumeWeight.C new file mode 100644 index 0000000000000000000000000000000000000000..9cd8e417f1c497ee21f638702c6d7b56edcd3553 --- /dev/null +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcCellVolumeWeight.C @@ -0,0 +1,301 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H" +#include "tetOverlapVolume.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::meshToMeshNew::calcCellVolumeWeight +( + const polyMesh& src, + const polyMesh& tgt, + const label srcSeedI, + const label tgtSeedI, + const labelList& srcCellIDs, + boolList& mapFlag, + label& startSeedI +) +{ + label srcCellI = srcSeedI; + label tgtCellI = tgtSeedI; + + List<DynamicList<label> > srcToTgtAddr(src.nCells()); + List<DynamicList<scalar> > srcToTgtWght(src.nCells()); + + List<DynamicList<label> > tgtToSrcAddr(tgt.nCells()); + List<DynamicList<scalar> > tgtToSrcWght(tgt.nCells()); + + // list of tgt cell neighbour cells + DynamicList<label> nbrTgtCells(10); + + // list of tgt cells currently visited for srcCellI to avoid multiple hits + DynamicList<label> visitedTgtCells(10); + + // list to keep track of tgt cells used to seed src cells + labelList seedCells(src.nCells(), -1); + seedCells[srcCellI] = tgtCellI; + + const scalarField& srcVol = src.cellVolumes(); + + do + { + nbrTgtCells.clear(); + visitedTgtCells.clear(); + + // append initial target cell and neighbours + nbrTgtCells.append(tgtCellI); + appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells); + + do + { + tgtCellI = nbrTgtCells.remove(); + visitedTgtCells.append(tgtCellI); + + scalar vol = interVol(src, tgt, srcCellI, tgtCellI); + + // accumulate addressing and weights for valid intersection + if (vol/srcVol[srcCellI] > tolerance_) + { + // store src/tgt cell pair + srcToTgtAddr[srcCellI].append(tgtCellI); + srcToTgtWght[srcCellI].append(vol); + + tgtToSrcAddr[tgtCellI].append(srcCellI); + tgtToSrcWght[tgtCellI].append(vol); + + appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells); + + // accumulate intersection volume + V_ += vol; + } + } + while (!nbrTgtCells.empty()); + + mapFlag[srcCellI] = false; + + // find new source seed cell + setNextCells + ( + startSeedI, + srcCellI, + tgtCellI, + src, + tgt, + srcCellIDs, + mapFlag, + visitedTgtCells, + seedCells + ); + } + while (srcCellI != -1); + + // transfer addressing into persistent storage + forAll(srcToTgtCellAddr_, i) + { + srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]); + srcToTgtCellWght_[i].transfer(srcToTgtWght[i]); + } + + forAll(tgtToSrcCellAddr_, i) + { + tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]); + tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]); + } +} + + +void Foam::meshToMeshNew::setNextCells +( + label& startSeedI, + label& srcCellI, + label& tgtCellI, + const polyMesh& src, + const polyMesh& tgt, + const labelList& srcCellIDs, + const boolList& mapFlag, + const DynamicList<label>& visitedCells, + labelList& seedCells +) const +{ + const labelList& srcNbrCells = src.cellCells()[srcCellI]; + + // set possible seeds for later use by querying all src cell neighbours + // with all visited target cells + bool valuesSet = false; + forAll(srcNbrCells, i) + { + label cellS = srcNbrCells[i]; + + if (mapFlag[cellS] && seedCells[cellS] == -1) + { + forAll(visitedCells, j) + { + label cellT = visitedCells[j]; + + if (intersect(src, tgt, cellS, cellT)) + { + seedCells[cellS] = cellT; + + if (!valuesSet) + { + srcCellI = cellS; + tgtCellI = cellT; + valuesSet = true; + } + } + } + } + } + + // set next src and tgt cells if not set above + if (valuesSet) + { + return; + } + else + { + // try to use existing seed + bool foundNextSeed = false; + for (label i = startSeedI; i < srcCellIDs.size(); i++) + { + label cellS = srcCellIDs[i]; + + if (mapFlag[cellS]) + { + if (!foundNextSeed) + { + startSeedI = i; + foundNextSeed = true; + } + + if (seedCells[cellS] != -1) + { + srcCellI = cellS; + tgtCellI = seedCells[cellS]; + + return; + } + } + } + + // perform new search to find match + if (debug) + { + Pout<< "Advancing front stalled: searching for new " + << "target cell" << endl; + } + + bool restart = + findInitialSeeds + ( + src, + tgt, + srcCellIDs, + mapFlag, + startSeedI, + srcCellI, + tgtCellI + ); + + if (restart) + { + // successfully found new starting seed-pair + return; + } + } + + // if we have got to here, there are no more src/tgt cell intersections + srcCellI = -1; + tgtCellI = -1; +} + + +bool Foam::meshToMeshNew::intersect +( + const polyMesh& src, + const polyMesh& tgt, + const label srcCellI, + const label tgtCellI +) const +{ + scalar threshold = tolerance_*src.cellVolumes()[srcCellI]; + + tetOverlapVolume overlapEngine; + + treeBoundBox bbTgtCell + ( + pointField + ( + tgt.points(), + tgt.cellPoints()[tgtCellI] + ) + ); + + return overlapEngine.cellCellOverlapMinDecomp + ( + src, + srcCellI, + tgt, + tgtCellI, + bbTgtCell, + threshold + ); +} + + +Foam::scalar Foam::meshToMeshNew::interVol +( + const polyMesh& src, + const polyMesh& tgt, + const label srcCellI, + const label tgtCellI +) const +{ + tetOverlapVolume overlapEngine; + + treeBoundBox bbTgtCell + ( + pointField + ( + tgt.points(), + tgt.cellPoints()[tgtCellI] + ) + ); + + scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp + ( + src, + srcCellI, + tgt, + tgtCellI, + bbTgtCell + ); + + return vol; +} + + +// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C new file mode 100644 index 0000000000000000000000000000000000000000..8e6a17511ef590af5f8eddfa339a73994a146109 --- /dev/null +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C @@ -0,0 +1,159 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::meshToMeshNew::calcDirect +( + const polyMesh& src, + const polyMesh& tgt, + const label srcSeedI, + const label tgtSeedI +) +{ + // store a list of src cells already mapped + boolList srcSeedFlag(src.nCells(), true); + labelList srcTgtSeed(src.nCells(), -1); + + List<DynamicList<label> > srcToTgt(src.nCells()); + List<DynamicList<label> > tgtToSrc(tgt.nCells()); + + DynamicList<label> srcSeeds; + + const scalarField& srcVc = src.cellVolumes(); + const scalarField& tgtVc = tgt.cellVolumes(); + + label srcCellI = srcSeedI; + label tgtCellI = tgtSeedI; + + do + { + // store src/tgt cell pair + srcToTgt[srcCellI].append(tgtCellI); + tgtToSrc[tgtCellI].append(srcCellI); + + // mark source cell srcSeedI as matched + srcSeedFlag[srcCellI] = false; + + // accumulate intersection volume + V_ += srcVc[srcCellI]; + + // find new source seed cell + appendToDirectSeeds + ( + src, + tgt, + srcSeedFlag, + srcTgtSeed, + srcSeeds, + srcCellI, + tgtCellI + ); + } + while (srcCellI >= 0); + + // transfer addressing into persistent storage + forAll(srcToTgtCellAddr_, i) + { + scalar v = srcVc[i]; + srcToTgtCellAddr_[i].transfer(srcToTgt[i]); + srcToTgtCellWght_[i] = scalarList(srcToTgtCellAddr_[i].size(), v); + } + + forAll(tgtToSrcCellAddr_, i) + { + scalar v = tgtVc[i]; + tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]); + tgtToSrcCellWght_[i] = scalarList(tgtToSrcCellAddr_[i].size(), v); + } +} + + +void Foam::meshToMeshNew::appendToDirectSeeds +( + const polyMesh& src, + const polyMesh& tgt, + boolList& mapFlag, + labelList& srcTgtSeed, + DynamicList<label>& srcSeeds, + label& srcSeedI, + label& tgtSeedI +) const +{ + const labelList& srcNbr = src.cellCells()[srcSeedI]; + const labelList& tgtNbr = tgt.cellCells()[tgtSeedI]; + + const vectorField& srcCentre = src.cellCentres(); + + forAll(srcNbr, i) + { + label srcI = srcNbr[i]; + + if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1)) + { + // source cell srcI not yet mapped + + // identfy if target cell exists for source cell srcI + bool found = false; + forAll(tgtNbr, j) + { + label tgtI = tgtNbr[j]; + + if (tgt.pointInCell(srcCentre[srcI], tgtI)) + { + // new match - append to lists + found = true; + + srcTgtSeed[srcI] = tgtI; + srcSeeds.append(srcI); + + break; + } + } + + if (!found) + { + // no match available for source cell srcI + mapFlag[srcI] = false; + } + } + } + + if (srcSeeds.size()) + { + srcSeedI = srcSeeds.remove(); + tgtSeedI = srcTgtSeed[srcSeedI]; + } + else + { + srcSeedI = -1; + tgtSeedI = -1; + } +} + + +// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMapNearest.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMapNearest.C new file mode 100644 index 0000000000000000000000000000000000000000..0276ee560495f8679eda7c7aa7b7710a41030aa3 --- /dev/null +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMapNearest.C @@ -0,0 +1,265 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012-2013 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 "meshToMeshNew.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::meshToMeshNew::calcMapNearest +( + const polyMesh& src, + const polyMesh& tgt, + const label srcSeedI, + const label tgtSeedI, + const labelList& srcCellIDs, + boolList& mapFlag, + label& startSeedI +) +{ + List<DynamicList<label> > srcToTgt(src.nCells()); + List<DynamicList<label> > tgtToSrc(tgt.nCells()); + + const scalarField& srcVc = src.cellVolumes(); + const scalarField& tgtVc = tgt.cellVolumes(); + + label srcCellI = srcSeedI; + label tgtCellI = tgtSeedI; + + do + { + // find nearest tgt cell + findNearestCell(src, tgt, srcCellI, tgtCellI); + + // store src/tgt cell pair + srcToTgt[srcCellI].append(tgtCellI); + tgtToSrc[tgtCellI].append(srcCellI); + + // mark source cell srcCellI and tgtCellI as matched + mapFlag[srcCellI] = false; + + // accumulate intersection volume + V_ += srcVc[srcCellI]; + + // find new source cell + setNextNearestCells + ( + startSeedI, + srcCellI, + tgtCellI, + mapFlag, + src, + tgt, + srcCellIDs + ); + } + while (srcCellI >= 0); + + + // for the case of multiple source cells per target cell, select the + // nearest source cell only and discard the others + const vectorField& srcCc = src.cellCentres(); + const vectorField& tgtCc = tgt.cellCentres(); + + forAll(tgtToSrc, targetCellI) + { + if (tgtToSrc[targetCellI].size() > 1) + { + const vector& tgtC = tgtCc[tgtCellI]; + + DynamicList<label>& srcCells = tgtToSrc[targetCellI]; + + label srcCellI = srcCells[0]; + scalar d = magSqr(tgtC - srcCc[srcCellI]); + + for (label i = 1; i < srcCells.size(); i++) + { + label srcI = srcCells[i]; + scalar dNew = magSqr(tgtC - srcCc[srcI]); + if (dNew < d) + { + d = dNew; + srcCellI = srcI; + } + } + + srcCells.clear(); + srcCells.append(srcCellI); + } + } + + // If there are more target cells than source cells, some target cells + // might not yet be mapped + forAll(tgtToSrc, tgtCellI) + { + if (tgtToSrc[tgtCellI].empty()) + { + label srcCellI = findMappedSrcCell(tgt, tgtCellI, tgtToSrc); + + findNearestCell(tgt, src, tgtCellI, srcCellI); + + tgtToSrc[tgtCellI].append(srcCellI); + } + } + + // transfer addressing into persistent storage + forAll(srcToTgtCellAddr_, i) + { + scalar v = srcVc[i]; + srcToTgtCellWght_[i] = scalarList(srcToTgt[i].size(), v); + srcToTgtCellAddr_[i].transfer(srcToTgt[i]); + } + + forAll(tgtToSrcCellAddr_, i) + { + scalar v = tgtVc[i]; + tgtToSrcCellWght_[i] = scalarList(tgtToSrc[i].size(), v); + tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]); + } +} + + +void Foam::meshToMeshNew::findNearestCell +( + const polyMesh& src, + const polyMesh& tgt, + const label srcCellI, + label& tgtCellI +) +{ + const vectorField& srcC = src.cellCentres(); + const vectorField& tgtC = tgt.cellCentres(); + + const vector& srcP = srcC[srcCellI]; + + DynamicList<label> tgtCells(10); + tgtCells.append(tgtCellI); + + DynamicList<label> visitedCells(10); + + scalar d = GREAT; + + do + { + label tgtI = tgtCells.remove(); + visitedCells.append(tgtI); + + scalar dTest = magSqr(tgtC[tgtI] - srcP); + if (dTest < d) + { + tgtCellI = tgtI; + d = dTest; + appendNbrCells(tgtCellI, tgt, visitedCells, tgtCells); + } + + } while (tgtCells.size() > 0); +} + + +void Foam::meshToMeshNew::setNextNearestCells +( + label& startSeedI, + label& srcCellI, + label& tgtCellI, + boolList& mapFlag, + const polyMesh& src, + const polyMesh& tgt, + const labelList& srcCellIDs +) +{ + const labelList& srcNbr = src.cellCells()[srcCellI]; + + srcCellI = -1; + forAll(srcNbr, i) + { + label cellI = srcNbr[i]; + if (mapFlag[cellI]) + { + srcCellI = cellI; + startSeedI = cellI + 1; + + return; + } + } + + (void)findInitialSeeds + ( + src, + tgt, + srcCellIDs, + mapFlag, + startSeedI, + srcCellI, + tgtCellI + ); +} + + +Foam::label Foam::meshToMeshNew::findMappedSrcCell +( + const polyMesh& tgt, + const label tgtCellI, + const List<DynamicList<label> >& tgtToSrc +) const +{ + DynamicList<label> testCells(10); + DynamicList<label> visitedCells(10); + + testCells.append(tgtCellI); + + do + { + // search target tgtCellI neighbours for match with source cell + label tgtI = testCells.remove(); + + if (findIndex(visitedCells, tgtI) == -1) + { + visitedCells.append(tgtI); + + if (tgtToSrc[tgtI].size()) + { + return tgtToSrc[tgtI][0]; + } + else + { + const labelList& nbrCells = tgt.cellCells()[tgtI]; + + forAll(nbrCells, i) + { + if (findIndex(visitedCells, nbrCells[i]) == -1) + { + testCells.append(nbrCells[i]); + } + } + } + } + } while (testCells.size()); + + // did not find any match - should not be possible to get here! + return -1; +} + + +// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C index 9ecda6cb4fb537e1cd349501d775e055edccad30..4add751668f357bbc37cbd493660df5ba5820500 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C @@ -29,7 +29,6 @@ License #include "globalIndex.H" #include "mergePoints.H" #include "treeBoundBox.H" -#include "tetOverlapVolume.H" #include "indexedOctree.H" #include "treeDataCell.H" #include "ListOps.H" @@ -44,14 +43,15 @@ namespace Foam const char* Foam::NamedEnum < Foam::meshToMeshNew::interpolationMethod, - 2 + 3 >::names[] = { - "map", + "direct", + "mapNearest", "cellVolumeWeight" }; - const NamedEnum<meshToMeshNew::interpolationMethod, 2> + const NamedEnum<meshToMeshNew::interpolationMethod, 3> meshToMeshNew::interpolationMethodNames_; } @@ -190,130 +190,30 @@ bool Foam::meshToMeshNew::findInitialSeeds } -void Foam::meshToMeshNew::appendToDirectSeeds +void Foam::meshToMeshNew::appendNbrCells ( - const polyMesh& src, - const polyMesh& tgt, - boolList& mapFlag, - labelList& srcTgtSeed, - DynamicList<label>& srcSeeds, - label& srcSeedI, - label& tgtSeedI + const label cellI, + const polyMesh& mesh, + const DynamicList<label>& visitedCells, + DynamicList<label>& nbrCellIDs ) const { - const labelList& srcNbr = src.cellCells()[srcSeedI]; - const labelList& tgtNbr = tgt.cellCells()[tgtSeedI]; - - const vectorField& srcCentre = src.cellCentres(); + const labelList& nbrCells = mesh.cellCells()[cellI]; - forAll(srcNbr, i) + // filter out cells already visited from cell neighbours + forAll(nbrCells, i) { - label srcI = srcNbr[i]; + label nbrCellI = nbrCells[i]; - if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1)) + if + ( + (findIndex(visitedCells, nbrCellI) == -1) + && (findIndex(nbrCellIDs, nbrCellI) == -1) + ) { - // source cell srcI not yet mapped - - // identfy if target cell exists for source cell srcI - bool found = false; - forAll(tgtNbr, j) - { - label tgtI = tgtNbr[j]; - - if (tgt.pointInCell(srcCentre[srcI], tgtI)) - { - // new match - append to lists - found = true; - - srcTgtSeed[srcI] = tgtI; - srcSeeds.append(srcI); - - break; - } - } - - if (!found) - { - // no match available for source cell srcI - mapFlag[srcI] = false; - } + nbrCellIDs.append(nbrCellI); } } - - if (srcSeeds.size()) - { - srcSeedI = srcSeeds.remove(); - tgtSeedI = srcTgtSeed[srcSeedI]; - } - else - { - srcSeedI = -1; - tgtSeedI = -1; - } -} - - -void Foam::meshToMeshNew::calcDirect -( - const polyMesh& src, - const polyMesh& tgt, - const label srcSeedI, - const label tgtSeedI -) -{ - // store a list of src cells already mapped - boolList srcSeedFlag(src.nCells(), true); - labelList srcTgtSeed(src.nCells(), -1); - - List<DynamicList<label> > srcToTgt(src.nCells()); - List<DynamicList<label> > tgtToSrc(tgt.nCells()); - - DynamicList<label> srcSeeds; - - const scalarField& srcVc = src.cellVolumes(); - - label srcCellI = srcSeedI; - label tgtCellI = tgtSeedI; - - do - { - // store src/tgt cell pair - srcToTgt[srcCellI].append(tgtCellI); - tgtToSrc[tgtCellI].append(srcCellI); - - // mark source cell srcSeedI as matched - srcSeedFlag[srcCellI] = false; - - // accumulate intersection volume - V_ += srcVc[srcCellI]; - - // find new source seed cell - appendToDirectSeeds - ( - src, - tgt, - srcSeedFlag, - srcTgtSeed, - srcSeeds, - srcCellI, - tgtCellI - ); - - } - while (srcCellI >= 0); - - // transfer addressing into persistent storage - forAll(srcToTgtCellAddr_, i) - { - srcToTgtCellAddr_[i].transfer(srcToTgt[i]); - srcToTgtCellWght_[i] = scalarList(srcToTgtCellAddr_[i].size(), 1.0); - } - - forAll(tgtToSrcCellAddr_, i) - { - tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]); - tgtToSrcCellWght_[i] = scalarList(tgtToSrcCellAddr_[i].size(), 1.0); - } } @@ -347,310 +247,13 @@ void Foam::meshToMeshNew::normaliseWeights maxW = max(maxW, s/Vc); } - Info<< type() << ": " << descriptor << " weights min/max = " + Info<< " " << descriptor << " weights min/max = " << returnReduce(minW, minOp<scalar>()) << ", " << returnReduce(maxW, maxOp<scalar>()) << endl; } } -void Foam::meshToMeshNew::appendNbrTgtCells -( - const label tgtCellI, - const polyMesh& tgt, - const DynamicList<label>& visitedTgtCells, - DynamicList<label>& nbrTgtCellIDs -) const -{ - const labelList& nbrCells = tgt.cellCells()[tgtCellI]; - - // filter out cells already visited from cell neighbours - forAll(nbrCells, i) - { - label nbrCellI = nbrCells[i]; - - if - ( - (findIndex(visitedTgtCells, nbrCellI) == -1) - && (findIndex(nbrTgtCellIDs, nbrCellI) == -1) - ) - { - nbrTgtCellIDs.append(nbrCellI); - } - } -} - - -void Foam::meshToMeshNew::setNextCells -( - label& startSeedI, - label& srcCellI, - label& tgtCellI, - const polyMesh& src, - const polyMesh& tgt, - const labelList& srcCellIDs, - const boolList& mapFlag, - const DynamicList<label>& visitedCells, - labelList& seedCells -) const -{ - const labelList& srcNbrCells = src.cellCells()[srcCellI]; - - // set possible seeds for later use by querying all src cell neighbours - // with all visited target cells - bool valuesSet = false; - forAll(srcNbrCells, i) - { - label cellS = srcNbrCells[i]; - - if (mapFlag[cellS] && seedCells[cellS] == -1) - { - forAll(visitedCells, j) - { - label cellT = visitedCells[j]; - - if (intersect(src, tgt, cellS, cellT)) - { - seedCells[cellS] = cellT; - - if (!valuesSet) - { - srcCellI = cellS; - tgtCellI = cellT; - valuesSet = true; - } - } - } - } - } - - // set next src and tgt cells if not set above - if (valuesSet) - { - return; - } - else - { - // try to use existing seed - bool foundNextSeed = false; - for (label i = startSeedI; i < srcCellIDs.size(); i++) - { - label cellS = srcCellIDs[i]; - - if (mapFlag[cellS]) - { - if (!foundNextSeed) - { - startSeedI = i; - foundNextSeed = true; - } - - if (seedCells[cellS] != -1) - { - srcCellI = cellS; - tgtCellI = seedCells[cellS]; - - return; - } - } - } - - // perform new search to find match - if (debug) - { - Pout<< "Advancing front stalled: searching for new " - << "target cell" << endl; - } - - bool restart = - findInitialSeeds - ( - src, - tgt, - srcCellIDs, - mapFlag, - startSeedI, - srcCellI, - tgtCellI - ); - - if (restart) - { - // successfully found new starting seed-pair - return; - } - } - - // if we have got to here, there are no more src/tgt cell intersections - srcCellI = -1; - tgtCellI = -1; -} - - -bool Foam::meshToMeshNew::intersect -( - const polyMesh& src, - const polyMesh& tgt, - const label srcCellI, - const label tgtCellI -) const -{ - scalar threshold = tolerance_*src.cellVolumes()[srcCellI]; - - tetOverlapVolume overlapEngine; - - treeBoundBox bbTgtCell - ( - pointField - ( - tgt.points(), - tgt.cellPoints()[tgtCellI] - ) - ); - - return overlapEngine.cellCellOverlapMinDecomp - ( - src, - srcCellI, - tgt, - tgtCellI, - bbTgtCell, - threshold - ); -} - - -Foam::scalar Foam::meshToMeshNew::interVol -( - const polyMesh& src, - const polyMesh& tgt, - const label srcCellI, - const label tgtCellI -) const -{ - tetOverlapVolume overlapEngine; - - treeBoundBox bbTgtCell - ( - pointField - ( - tgt.points(), - tgt.cellPoints()[tgtCellI] - ) - ); - - scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp - ( - src, - srcCellI, - tgt, - tgtCellI, - bbTgtCell - ); - - return vol; -} - - -void Foam::meshToMeshNew::calcIndirect -( - const polyMesh& src, - const polyMesh& tgt, - const label srcSeedI, - const label tgtSeedI, - const labelList& srcCellIDs, - boolList& mapFlag, - label& startSeedI -) -{ - label srcCellI = srcSeedI; - label tgtCellI = tgtSeedI; - - List<DynamicList<label> > srcToTgtAddr(src.nCells()); - List<DynamicList<scalar> > srcToTgtWght(src.nCells()); - - List<DynamicList<label> > tgtToSrcAddr(tgt.nCells()); - List<DynamicList<scalar> > tgtToSrcWght(tgt.nCells()); - - // list of tgt cell neighbour cells - DynamicList<label> nbrTgtCells(10); - - // list of tgt cells currently visited for srcCellI to avoid multiple hits - DynamicList<label> visitedTgtCells(10); - - // list to keep track of tgt cells used to seed src cells - labelList seedCells(src.nCells(), -1); - seedCells[srcCellI] = tgtCellI; - - const scalarField& srcVol = src.cellVolumes(); - - do - { - nbrTgtCells.clear(); - visitedTgtCells.clear(); - - // append initial target cell and neighbours - nbrTgtCells.append(tgtCellI); - appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells); - - do - { - tgtCellI = nbrTgtCells.remove(); - visitedTgtCells.append(tgtCellI); - - scalar vol = interVol(src, tgt, srcCellI, tgtCellI); - - // accumulate addressing and weights for valid intersection - if (vol/srcVol[srcCellI] > tolerance_) - { - // store src/tgt cell pair - srcToTgtAddr[srcCellI].append(tgtCellI); - srcToTgtWght[srcCellI].append(vol); - - tgtToSrcAddr[tgtCellI].append(srcCellI); - tgtToSrcWght[tgtCellI].append(vol); - - appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells); - - // accumulate intersection volume - V_ += vol; - } - } - while (!nbrTgtCells.empty()); - - mapFlag[srcCellI] = false; - - // find new source seed cell - setNextCells - ( - startSeedI, - srcCellI, - tgtCellI, - src, - tgt, - srcCellIDs, - mapFlag, - visitedTgtCells, - seedCells - ); - } - while (srcCellI != -1); - - // transfer addressing into persistent storage - forAll(srcToTgtCellAddr_, i) - { - srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]); - srcToTgtCellWght_[i].transfer(srcToTgtWght[i]); - } - - forAll(tgtToSrcCellAddr_, i) - { - tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]); - tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]); - } -} - - void Foam::meshToMeshNew::calcAddressing ( const polyMesh& src, @@ -723,14 +326,28 @@ void Foam::meshToMeshNew::calcAddressing switch (method_) { - case imMap: + case imDirect: { calcDirect(src, tgt, srcSeedI, tgtSeedI); break; } + case imMapNearest: + { + calcMapNearest + ( + src, + tgt, + srcSeedI, + tgtSeedI, + srcCellIDs, + mapFlag, + startSeedI + ); + break; + } case imCellVolumeWeight: { - calcIndirect + calcCellVolumeWeight ( src, tgt, @@ -765,42 +382,24 @@ void Foam::meshToMeshNew::calcAddressing } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::meshToMeshNew::meshToMeshNew -( - const polyMesh& src, - const polyMesh& tgt, - const interpolationMethod& method -) -: - srcRegionName_(src.name()), - tgtRegionName_(tgt.name()), - srcToTgtCellAddr_(), - tgtToSrcCellAddr_(), - srcToTgtCellWght_(), - tgtToSrcCellWght_(), - method_(method), - V_(0.0), - singleMeshProc_(-1), - srcMapPtr_(NULL), - tgtMapPtr_(NULL) +void Foam::meshToMeshNew::calculate() { - Info<< "Creating mesh-to-mesh addressing for " << src.name() - << " and " << tgt.name() << " regions" << endl; + Info<< "Creating mesh-to-mesh addressing for " << srcRegion_.name() + << " and " << tgtRegion_.name() << " regions using " + << interpolationMethodNames_[method_] << endl; - singleMeshProc_ = calcDistribution(src, tgt); + singleMeshProc_ = calcDistribution(srcRegion_, tgtRegion_); if (singleMeshProc_ == -1) { // create global indexing for src and tgt meshes - globalIndex globalSrcCells(src.nCells()); - globalIndex globalTgtCells(tgt.nCells()); + globalIndex globalSrcCells(srcRegion_.nCells()); + globalIndex globalTgtCells(tgtRegion_.nCells()); // Create processor map of overlapping cells. This map gets // (possibly remote) cells from the tgt mesh such that they (together) // cover all of the src mesh - autoPtr<mapDistribute> mapPtr = calcProcMap(src, tgt); + autoPtr<mapDistribute> mapPtr = calcProcMap(srcRegion_, tgtRegion_); const mapDistribute& map = mapPtr(); pointField newTgtPoints; @@ -812,7 +411,7 @@ Foam::meshToMeshNew::meshToMeshNew distributeAndMergeCells ( map, - tgt, + tgtRegion_, globalTgtCells, newTgtPoints, newTgtFaces, @@ -827,9 +426,9 @@ Foam::meshToMeshNew::meshToMeshNew ( IOobject ( - "newTgt::" + Foam::name(Pstream::myProcNo()), - tgt.time().timeName(), - tgt.time(), + "newTgt." + Foam::name(Pstream::myProcNo()), + tgtRegion_.time().timeName(), + tgtRegion_.time(), IOobject::NO_READ ), xferMove(newTgtPoints), @@ -862,9 +461,9 @@ Foam::meshToMeshNew::meshToMeshNew if (debug) { Pout<< "Created newTgt mesh:" << nl - << " old cells = " << tgt.nCells() + << " old cells = " << tgtRegion_.nCells() << ", new cells = " << newTgt.nCells() << nl - << " old faces = " << tgt.nFaces() + << " old faces = " << tgtRegion_.nFaces() << ", new faces = " << newTgt.nFaces() << endl; if (debug > 1) @@ -874,7 +473,7 @@ Foam::meshToMeshNew::meshToMeshNew } } - calcAddressing(src, newTgt); + calcAddressing(srcRegion_, newTgt); // per source cell the target cell address in newTgt mesh forAll(srcToTgtCellAddr_, i) @@ -901,7 +500,7 @@ Foam::meshToMeshNew::meshToMeshNew ( Pstream::nonBlocking, List<labelPair>(), - tgt.nCells(), + tgtRegion_.nCells(), map.constructMap(), map.subMap(), tgtToSrcCellAddr_, @@ -914,7 +513,7 @@ Foam::meshToMeshNew::meshToMeshNew ( Pstream::nonBlocking, List<labelPair>(), - tgt.nCells(), + tgtRegion_.nCells(), map.constructMap(), map.subMap(), tgtToSrcCellWght_, @@ -926,7 +525,7 @@ Foam::meshToMeshNew::meshToMeshNew normaliseWeights ( "source", - src.cellVolumes(), + srcRegion_.cellVolumes(), srcToTgtCellAddr_, srcToTgtCellWght_ ); @@ -934,7 +533,7 @@ Foam::meshToMeshNew::meshToMeshNew normaliseWeights ( "target", - tgt.cellVolumes(), + tgtRegion_.cellVolumes(), tgtToSrcCellAddr_, tgtToSrcCellWght_ ); @@ -955,12 +554,12 @@ Foam::meshToMeshNew::meshToMeshNew } else { - calcAddressing(src, tgt); + calcAddressing(srcRegion_, tgtRegion_); normaliseWeights ( "source", - src.cellVolumes(), + srcRegion_.cellVolumes(), srcToTgtCellAddr_, srcToTgtCellWght_ ); @@ -968,7 +567,7 @@ Foam::meshToMeshNew::meshToMeshNew normaliseWeights ( "target", - tgt.cellVolumes(), + tgtRegion_.cellVolumes(), tgtToSrcCellAddr_, tgtToSrcCellWght_ ); @@ -978,6 +577,162 @@ Foam::meshToMeshNew::meshToMeshNew } +const Foam::PtrList<Foam::AMIPatchToPatchInterpolation>& +Foam::meshToMeshNew::patchAMIs() const +{ + if (patchAMIs_.empty()) + { + 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() << endl; + + Info<< incrIndent; + + patchAMIs_.set + ( + i, + new AMIPatchToPatchInterpolation + ( + srcPP, + tgtPP, + faceAreaIntersect::tmMesh, + true // flip target patch since patch normals are aligned + ) + ); + + Info<< decrIndent; + } + } + + return patchAMIs_; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::meshToMeshNew::meshToMeshNew +( + const polyMesh& src, + const polyMesh& tgt, + const interpolationMethod& method, + bool interpAllPatches +) +: + srcRegion_(src), + tgtRegion_(tgt), + srcPatchID_(), + tgtPatchID_(), + patchAMIs_(), + 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(); + + if (srcBM.size() != tgtBM.size()) + { + FatalErrorIn + ( + "Foam::meshToMeshNew::meshToMeshNew" + "(" + "const polyMesh&, " + "const polyMesh&, " + "const interpolationMethod&" + ")" + ) << "Source and target meshes are dissimiar:" << nl + << " Source patches: " << srcBM.size() << nl + << " Target patches: " << tgtBM.size() << exit(FatalError); + } + + DynamicList<label> patchID(src.boundaryMesh().size()); + + forAll(srcBM, patchI) + { + const polyPatch& pp = srcBM[patchI]; + if (!polyPatch::constraintType(pp.type())) + { + patchID.append(pp.index()); + } + } + + srcPatchID_.transfer(patchID); + tgtPatchID_ = srcPatchID_; + } + + // calculate volume addressing and weights + calculate(); + + // calculate patch addressing and weights + (void)patchAMIs(); +} + + +Foam::meshToMeshNew::meshToMeshNew +( + const polyMesh& src, + const polyMesh& tgt, + const interpolationMethod& method, + const HashTable<word>& patchMap +) +: + srcRegion_(src), + tgtRegion_(tgt), + srcPatchID_(), + tgtPatchID_(), + patchAMIs_(), + srcToTgtCellAddr_(), + tgtToSrcCellAddr_(), + srcToTgtCellWght_(), + tgtToSrcCellWght_(), + method_(method), + V_(0.0), + singleMeshProc_(-1), + srcMapPtr_(NULL), + tgtMapPtr_(NULL) +{ + srcPatchID_.setSize(patchMap.size()); + tgtPatchID_.setSize(patchMap.size()); + + label i = 0; + forAllConstIter(HashTable<word>, patchMap, iter) + { + const word& srcPatchName = iter.key(); + const word& tgtPatchName = iter(); + + const polyPatch& srcPatch = srcRegion_.boundaryMesh()[srcPatchName]; + const polyPatch& tgtPatch = tgtRegion_.boundaryMesh()[tgtPatchName]; + + srcPatchID_[i] = srcPatch.index(); + tgtPatchID_[i] = tgtPatch.index(); + i++; + } + + // calculate volume addressing and weights + calculate(); + + // calculate patch addressing and weights + (void)patchAMIs(); +} + + // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // Foam::meshToMeshNew::~meshToMeshNew() diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H index 85f0c1332e9e5857008d54553a6fe38a0b6ff7f8..8404ce7b9fa077ecc07ea02671d4673dfeb6981c 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H @@ -27,7 +27,19 @@ Class Description Class to calculate the cell-addressing between two overlapping meshes + Three methods are currently available: + - direct : 1-to-1 mapping between meshes + - mapNearest : assign nearest cell values without interpolation + - cellVolumeWeight : volume consistent mapping + + The \c direct and \c cellVolumeWeight options are volume conservative, + whereas mapNearest is non-conservative. + + SourceFiles + calcDirect.C + calcMapNearest.C + calcCellVolumeWeight.C meshToMeshNew.C meshToMeshNewTemplates.C @@ -41,6 +53,7 @@ SourceFiles #include "mapDistribute.H" #include "volFieldsFwd.H" #include "NamedEnum.H" +#include "AMIPatchToPatchInterpolation.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -60,22 +73,32 @@ public: //- Enumeration specifying required accuracy enum interpolationMethod { - imMap, + imDirect, + imMapNearest, imCellVolumeWeight }; - static const NamedEnum<interpolationMethod, 2> + static const NamedEnum<interpolationMethod, 3> interpolationMethodNames_; private: // Private data - //- Name of source mesh region - const word srcRegionName_; + //- Reference to the source mesh + const polyMesh& srcRegion_; + + //- Reference to the target mesh + const polyMesh& tgtRegion_; - //- Name of target mesh region - const word tgtRegionName_; + //- List of target patch IDs per source patch (local index) + List<label> srcPatchID_; + + //- List of source patch IDs per target patch (local index) + List<label> tgtPatchID_; + + //- List of AMIs between source and target patches + mutable PtrList<AMIPatchToPatchInterpolation> patchAMIs_; //- Source to target cell addressing labelListList srcToTgtCellAddr_; @@ -139,9 +162,36 @@ private: label& tgtSeedI ) const; + //- Append target cell neihgbour cells to cellIDs list + void appendNbrCells + ( + const label tgtCellI, + const polyMesh& tgt, + const DynamicList<label>& visitedTgtCells, + DynamicList<label>& nbrTgtCellIDs + ) const; + + //- Normalise the interpolation weights + void normaliseWeights + ( + const word& descriptor, + const scalarField& cellVolumes, + const labelListList& addr, + scalarListList& wght + ) const; + // Direct (one-to-one) mapping + //- Main driver routine for direct mapping + void calcDirect + ( + const polyMesh& src, + const polyMesh& tgt, + const label srcSeedI, + const label tgtSeedI + ); + //- Append to list of src mesh seed indices void appendToDirectSeeds ( @@ -154,35 +204,63 @@ private: label& tgtSeedI ) const; - //- Main driver routine for direct mapping - void calcDirect + // Nearest (non-conformal) mapping + + //- Main driver routine for nearest-mapping routine + void calcMapNearest ( const polyMesh& src, const polyMesh& tgt, const label srcSeedI, - const label tgtSeedI + const label tgtSeedI, + const labelList& srcCellIDs, + boolList& mapFlag, + label& startSeedI ); + //- Find target cell index of cell closest to source cell + void findNearestCell + ( + const polyMesh& src, + const polyMesh& tgt, + const label srcCellI, + label& tgtCellI + ); - // Indirect (non-conformal) mapping + //- Set the next pair of cells + void setNextNearestCells + ( + label& startSeedI, + label& srcCellI, + label& tgtCellI, + boolList& mapFlag, + const polyMesh& src, + const polyMesh& tgt, + const labelList& srcCellIDs + ); - //- Normalise the interpolation weights - void normaliseWeights + //- Find source cell for target cell + label findMappedSrcCell ( - const word& descriptor, - const scalarField& cellVolumes, - const labelListList& addr, - scalarListList& wght + const polyMesh& tgt, + const label tgtCellI, + const List<DynamicList<label> >& tgtToSrc ) const; - //- Append target cell neihgbour cells to cellIDs list - void appendNbrTgtCells + + // Cell volume weighted (non-conformal) interpolation + + //- Main driver routine for cell volume weighted interpolation + void calcCellVolumeWeight ( - const label tgtCellI, + const polyMesh& src, const polyMesh& tgt, - const DynamicList<label>& visitedTgtCells, - DynamicList<label>& nbrTgtCellIDs - ) const; + const label srcSeedI, + const label tgtSeedI, + const labelList& srcCellIDs, + boolList& mapFlag, + label& startSeedI + ); //- Set the next cells in the advancing front algorithm void setNextCells @@ -216,23 +294,17 @@ private: const label tgtCellI ) const; - //- Main driver routine for indirect mapping - void calcIndirect - ( - const polyMesh& src, - const polyMesh& tgt, - const label srcSeedI, - const label tgtSeedI, - const labelList& srcCellIDs, - boolList& mapFlag, - label& startSeedI - ); - //- Calculate the addressing between overalping regions of src and tgt - // meshes - main driver function + // meshes void calcAddressing(const polyMesh& src, const polyMesh& tgt); + //- Calculate - main driver function + void calculate(); + + //- Return the list of AMIs between source and target patches + const PtrList<AMIPatchToPatchInterpolation>& patchAMIs() const; + // Parallel operations @@ -306,7 +378,18 @@ public: ( const polyMesh& src, const polyMesh& tgt, - const interpolationMethod& method + const interpolationMethod& method, + const bool interpAllPatches = true + ); + + + //- Construct from source and target meshes + meshToMeshNew + ( + const polyMesh& src, + const polyMesh& tgt, + const interpolationMethod& method, + const HashTable<word>& patchMap ); @@ -318,6 +401,12 @@ public: // Access + //- Return const access to the source mesh + inline const polyMesh& srcRegion() const; + + //- Return const access to the target mesh + inline const polyMesh& tgtRegion() const; + //- Return const access to the source to target cell addressing inline const labelListList& srcToTgtCellAddr() const; @@ -432,60 +521,103 @@ public: ) const; - // Volume field mapping + // Source-to-target volume field mapping - //- Interpolare a field with a defined operation. Values + //- Interpolate a field with a defined operation. Values // passed in via 'result' are used to initialise the return - // value. Optionally interpolate patch values + // value template<class Type, class CombineOp> - void interpolate + void mapSrcToTgt ( const GeometricField<Type, fvPatchField, volMesh>& field, const CombineOp& cop, - GeometricField<Type, fvPatchField, volMesh>& result, - const bool interpPatches = true + GeometricField<Type, fvPatchField, volMesh>& result ) const; - //- Interpolare a field with a defined operation. The initial - // values of the result are set to zero. Optionally - // interpolate patch values + //- Interpolate a field with a defined operation. The initial + // values of the result are set to zero template<class Type, class CombineOp> - tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate + tmp<GeometricField<Type, fvPatchField, volMesh> > mapSrcToTgt ( const GeometricField<Type, fvPatchField, volMesh>& field, - const CombineOp& cop, - const bool interpPatches = true + const CombineOp& cop ) const; - //- Interpolare a tmp field with a defined operation. The - // initial values of the result are set to zero. Optionally - // interpolate patch values + //- Interpolate a tmp field with a defined operation. The + // initial values of the result are set to zero template<class Type, class CombineOp> - tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate + tmp<GeometricField<Type, fvPatchField, volMesh> > mapSrcToTgt ( const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield, - const CombineOp& cop, - const bool interpPatches = true + const CombineOp& cop ) const; //- Convenience function to map a field with a default - // operation (plusEqOp). Optionally interpolate patch values + // operation (plusEqOp) template<class Type> - tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate + tmp<GeometricField<Type, fvPatchField, volMesh> > mapSrcToTgt ( - const GeometricField<Type, fvPatchField, volMesh>& field, - const bool interpPatches = true + const GeometricField<Type, fvPatchField, volMesh>& field ) const; //- Convenience function to map a tmp field with a default - // operation (plusEqOp). Optionally interpolate patch values + // operation (plusEqOp) template<class Type> - tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate + tmp<GeometricField<Type, fvPatchField, volMesh> > mapSrcToTgt + ( + const tmp<GeometricField<Type, fvPatchField, volMesh> >& + tfield + ) const; + + + // Target-to-source volume field mapping + + //- Interpolate a field with a defined operation. Values + // passed in via 'result' are used to initialise the return + // value + template<class Type, class CombineOp> + void mapTgtToSrc + ( + const GeometricField<Type, fvPatchField, volMesh>& field, + const CombineOp& cop, + GeometricField<Type, fvPatchField, volMesh>& result + ) const; + + //- Interpolate a field with a defined operation. The initial + // values of the result are set to zero + template<class Type, class CombineOp> + tmp<GeometricField<Type, fvPatchField, volMesh> > mapTgtToSrc + ( + const GeometricField<Type, fvPatchField, volMesh>& field, + const CombineOp& cop + ) const; + + //- Interpolate a tmp field with a defined operation. The + // initial values of the result are set to zero + template<class Type, class CombineOp> + tmp<GeometricField<Type, fvPatchField, volMesh> > mapTgtToSrc ( const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield, - const bool interpPatches = true + const CombineOp& cop + ) const; + + //- Convenience function to map a field with a default + // operation (plusEqOp) + template<class Type> + tmp<GeometricField<Type, fvPatchField, volMesh> > mapTgtToSrc + ( + const GeometricField<Type, fvPatchField, volMesh>& field + ) const; + + //- Convenience function to map a tmp field with a default + // operation (plusEqOp) + template<class Type> + tmp<GeometricField<Type, fvPatchField, volMesh> > mapTgtToSrc + ( + const tmp<GeometricField<Type, fvPatchField, volMesh> >& + tfield ) const; }; diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H index 41068e8c455e4e9f5b8d7737acc2fde589df618c..157a9efd653b360d779425dd2cf21a0a855f00cf 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H @@ -27,6 +27,18 @@ License // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +inline const Foam::polyMesh& Foam::meshToMeshNew::srcRegion() const +{ + return srcRegion_; +} + + +inline const Foam::polyMesh& Foam::meshToMeshNew::tgtRegion() const +{ + return tgtRegion_; +} + + inline const Foam::labelListList& Foam::meshToMeshNew::srcToTgtCellAddr() const { diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C index 69bbc69667b466985a055d3c424b925eec0ca38f..dde5218fd485f78e852cde335dafb752eba8082a 100644 --- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C +++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C @@ -25,7 +25,6 @@ License #include "fvMesh.H" #include "volFields.H" -//#include "ops.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -56,31 +55,6 @@ namespace Foam } } }; - - //- Combine operator for maps/interpolations - template<class Type, class CombineOp> - class combineBinaryOp - { - const CombineOp& cop_; - - public: - - combineBinaryOp(const CombineOp& cop) - : - cop_(cop) - {} - - void operator() - ( - Type& x, - const label faceI, - const Type& y, - const scalar weight - ) const - { - cop_(x, weight*y); - } - }; } @@ -119,13 +93,13 @@ void Foam::meshToMeshNew::mapSrcToTgt "List<Type>&" ") const" ) << "Supplied field size is not equal to target mesh size" << nl - << " source mesh = " << srcToTgtCellAddr_.size() << nl - << " target mesh = " << tgtToSrcCellAddr_.size() << nl + << " source mesh = " << srcToTgtCellAddr_.size() << nl + << " target mesh = " << tgtToSrcCellAddr_.size() << nl << " supplied field = " << result.size() << abort(FatalError); } - combineBinaryOp<Type, CombineOp> cbop(cop); + multiplyWeightedOp<Type, CombineOp> cbop(cop); if (singleMeshProc_ == -1) { @@ -247,13 +221,13 @@ void Foam::meshToMeshNew::mapTgtToSrc "List<Type>&" ") const" ) << "Supplied field size is not equal to source mesh size" << nl - << " source mesh = " << srcToTgtCellAddr_.size() << nl - << " target mesh = " << tgtToSrcCellAddr_.size() << nl + << " source mesh = " << srcToTgtCellAddr_.size() << nl + << " target mesh = " << tgtToSrcCellAddr_.size() << nl << " supplied field = " << result.size() << abort(FatalError); } - combineBinaryOp<Type, CombineOp> cbop(cop); + multiplyWeightedOp<Type, CombineOp> cbop(cop); if (singleMeshProc_ == -1) { @@ -269,7 +243,6 @@ void Foam::meshToMeshNew::mapTgtToSrc if (tgtAddress.size()) { -// result[cellI] = pTraits<Type>::zero; result[cellI] *= (1.0 - sum(tgtWeight)); forAll(tgtAddress, i) { @@ -289,7 +262,6 @@ void Foam::meshToMeshNew::mapTgtToSrc if (tgtAddress.size()) { -// result[cellI] = pTraits<Type>::zero; result[cellI] *= (1.0 - sum(tgtWeight)); forAll(tgtAddress, i) { @@ -357,140 +329,181 @@ Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc template<class Type, class CombineOp> -void Foam::meshToMeshNew::interpolate +void Foam::meshToMeshNew::mapSrcToTgt ( const GeometricField<Type, fvPatchField, volMesh>& field, const CombineOp& cop, - GeometricField<Type, fvPatchField, volMesh>& result, - const bool interpPatches + GeometricField<Type, fvPatchField, volMesh>& result ) const { - const fvMesh& mesh = field.mesh(); + // clear any previously stored values - if (mesh.name() == srcRegionName_) - { - mapSrcToTgt(field, cop, result.internalField()); - } - else if (mesh.name() == tgtRegionName_) + + mapSrcToTgt(field, cop, result.internalField()); + + const PtrList<AMIPatchToPatchInterpolation>& AMIList = patchAMIs(); + + forAll(AMIList, i) { - mapTgtToSrc(field, cop, result.internalField()); + label srcPatchI = srcPatchID_[i]; + label tgtPatchI = tgtPatchID_[i]; + + const Field<Type>& srcField = field.boundaryField()[srcPatchI]; + Field<Type>& tgtField = result.boundaryField()[tgtPatchI]; + + tgtField = pTraits<Type>::zero; + + AMIList[i].interpolateToTarget + ( + srcField, + multiplyWeightedOp<Type, CombineOp>(cop), + tgtField + ); } - else - { - FatalErrorIn +} + + +template<class Type, class CombineOp> +Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > +Foam::meshToMeshNew::mapSrcToTgt +( + const GeometricField<Type, fvPatchField, volMesh>& field, + const CombineOp& cop +) const +{ + typedef GeometricField<Type, fvPatchField, volMesh> fieldType; + + const fvMesh& tgtMesh = static_cast<const fvMesh&>(tgtRegion_); + + tmp<fieldType> tresult + ( + new fieldType ( - "void Foam::meshToMeshNew::interpolate" - "(" - "const GeometricField<Type, fvPatchField, volMesh>&, " - "const CombineOp&, " - "GeometricField<Type, fvPatchField, volMesh>&, " - "const bool" - ") const" + IOobject + ( + type() + ".interpolate(" + field.name() + ")", + tgtMesh.time().timeName(), + tgtMesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + tgtMesh, + dimensioned<Type> + ( + "zero", + field.dimensions(), + pTraits<Type>::zero + ) ) - << "Supplied field " << field.name() << " did not originate from " - << "either the source or target meshes used to create this " - << "interpolation object" - << abort(FatalError); - } + ); - if (interpPatches) - { - switch (method_) - { - case imMap: - { - result.boundaryField() == field.boundaryField(); - break; - } - default: - { - notImplemented - ( - "void Foam::meshToMeshNew::interpolate" - "(" - "const GeometricField<Type, fvPatchField, volMesh>&, " - "const CombineOp&, " - "GeometricField<Type, fvPatchField, volMesh>&, " - "const bool" - ") const - non-conformal patches" - ) - - // do something... - } - } - } + mapSrcToTgt(field, cop, tresult()); + + return tresult; } template<class Type, class CombineOp> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > -Foam::meshToMeshNew::interpolate +Foam::meshToMeshNew::mapSrcToTgt +( + const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield, + const CombineOp& cop +) const +{ + return mapSrcToTgt(tfield(), cop); +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > +Foam::meshToMeshNew::mapSrcToTgt +( + const GeometricField<Type, fvPatchField, volMesh>& field +) const +{ + return mapSrcToTgt(field, plusEqOp<Type>()); +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > +Foam::meshToMeshNew::mapSrcToTgt +( + const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield +) const +{ + return mapSrcToTgt(tfield(), plusEqOp<Type>()); +} + + +template<class Type, class CombineOp> +void Foam::meshToMeshNew::mapTgtToSrc ( const GeometricField<Type, fvPatchField, volMesh>& field, const CombineOp& cop, - const bool interpPatches + GeometricField<Type, fvPatchField, volMesh>& result ) const { - typedef GeometricField<Type, fvPatchField, volMesh> fieldType; - - const fvMesh& mesh = field.mesh(); + mapTgtToSrc(field, cop, result.internalField()); - tmp<fieldType> tresult; + const PtrList<AMIPatchToPatchInterpolation>& AMIList = patchAMIs(); - if (mesh.name() == srcRegionName_) + forAll(AMIList, i) { - const fvMesh& tgtMesh = - mesh.time().lookupObject<fvMesh>(tgtRegionName_); + label srcPatchI = srcPatchID_[i]; + label tgtPatchI = tgtPatchID_[i]; - tresult = - new fieldType - ( - IOobject - ( - type() + "::interpolate(" + field.name() + ")", - tgtMesh.time().timeName(), - tgtMesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - tgtMesh, - dimensioned<Type> - ( - "zero", - field.dimensions(), - pTraits<Type>::zero - ) - ); - - interpolate(field, cop, tresult(), interpPatches); + Field<Type>& srcField = result.boundaryField()[srcPatchI]; + const Field<Type>& tgtField = field.boundaryField()[tgtPatchI]; + + srcField = pTraits<Type>::zero; + + AMIList[i].interpolateToSource + ( + tgtField, + multiplyWeightedOp<Type, CombineOp>(cop), + srcField + ); } - else if (mesh.name() == tgtRegionName_) - { - const fvMesh& srcMesh = - mesh.time().lookupObject<fvMesh>(srcRegionName_); +} + + +template<class Type, class CombineOp> +Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > +Foam::meshToMeshNew::mapTgtToSrc +( + const GeometricField<Type, fvPatchField, volMesh>& field, + const CombineOp& cop +) const +{ + typedef GeometricField<Type, fvPatchField, volMesh> fieldType; - tresult = - new fieldType + const fvMesh& srcMesh = static_cast<const fvMesh&>(srcRegion_); + + tmp<fieldType> tresult + ( + new fieldType + ( + IOobject ( - IOobject - ( - type() + "::interpolate(" + field.name() + ")", - srcMesh.time().timeName(), - srcMesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), + type() + ".interpolate(" + field.name() + ")", + srcMesh.time().timeName(), srcMesh, - dimensioned<Type> - ( - "zero", - field.dimensions(), - pTraits<Type>::zero - ) - ); - - interpolate(field, cop, tresult(), interpPatches); - } + IOobject::NO_READ, + IOobject::NO_WRITE + ), + srcMesh, + dimensioned<Type> + ( + "zero", + field.dimensions(), + pTraits<Type>::zero + ) + ) + ); + + mapTgtToSrc(field, cop, tresult()); return tresult; } @@ -498,44 +511,35 @@ Foam::meshToMeshNew::interpolate template<class Type, class CombineOp> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > -Foam::meshToMeshNew::interpolate +Foam::meshToMeshNew::mapTgtToSrc ( const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield, - const CombineOp& cop, - const bool interpPatches + const CombineOp& cop ) const { - return - interpolate - ( - tfield(), - combineBinaryOp<Type, CombineOp>(cop), - interpPatches - ); + return mapTgtToSrc(tfield(), cop); } template<class Type> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > -Foam::meshToMeshNew::interpolate +Foam::meshToMeshNew::mapTgtToSrc ( - const GeometricField<Type, fvPatchField, volMesh>& field, - const bool interpPatches + const GeometricField<Type, fvPatchField, volMesh>& field ) const { - return interpolate(field, plusEqOp<Type>(), interpPatches); + return mapTgtToSrc(field, plusEqOp<Type>()); } template<class Type> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> > -Foam::meshToMeshNew::interpolate +Foam::meshToMeshNew::mapTgtToSrc ( - const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield, - const bool interpPatches + const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield ) const { - return interpolate(tfield(), plusEqOp<Type>(), interpPatches); + return mapTgtToSrc(tfield(), plusEqOp<Type>()); } diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C index ba9416b6ed4a958e8c7b4705965464e48f5df1f4..90e6977a776db4ae4331a34eb6033c5ebae1c909 100644 --- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C @@ -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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -191,6 +191,35 @@ externalWallHeatFluxTemperatureFvPatchScalarField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::autoMap +( + const fvPatchFieldMapper& m +) +{ + mixedFvPatchScalarField::autoMap(m); + q_.autoMap(m); + h_.autoMap(m); + Ta_.autoMap(m); +} + + +void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::rmap +( + const fvPatchScalarField& ptf, + const labelList& addr +) +{ + mixedFvPatchScalarField::rmap(ptf, addr); + + const externalWallHeatFluxTemperatureFvPatchScalarField& tiptf = + refCast<const externalWallHeatFluxTemperatureFvPatchScalarField>(ptf); + + q_.rmap(tiptf.q_, addr); + h_.rmap(tiptf.h_, addr); + Ta_.rmap(tiptf.Ta_, addr); +} + + void Foam::externalWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs() { if (updated()) diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H index 44e96635bacdf6056c9b724dff868e8660977e9a..5fdaa446eecfa612c097f59a413cdbca817d4dca 100644 --- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -185,6 +185,22 @@ public: // Member functions + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchScalarField&, + const labelList& + ); + + // Evaluation functions //- Update the coefficients associated with the patch field diff --git a/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict b/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict index f70ad3e89bda2de93f8d64451856d211af71e2c9..1ee58bdb3e6049364e4dcefdc221f66701d298fb 100644 --- a/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict +++ b/tutorials/incompressible/pimpleFoam/TJunctionFan/system/createBafflesDict @@ -41,14 +41,16 @@ baffles name baffles; type wall; + XXX 9.8; patchFields { + epsilon { type epsilonWallFunction; Cmu 0.09; kappa 0.41; - E 9.8; + E $Cmu; // XXX; //9.8; value uniform 0; } k