From 6000d2fb97c50e7cc5828621c9010835fab57823 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 12 Apr 2013 15:52:12 +0100
Subject: [PATCH] ENH: Refectored mesh-to-mesh interpolation methods - now
run-time selectable
---
src/sampling/Make/files | 10 +-
.../meshToMeshNew/calcDirect.C | 159 -------
.../meshToMeshNew/calcMapNearest.C | 265 -----------
.../cellVolumeWeightMethod.C} | 241 ++++++----
.../cellVolumeWeight/cellVolumeWeightMethod.H | 138 ++++++
.../calcMethod/direct/directMethod.C | 303 +++++++++++++
.../calcMethod/direct/directMethod.H | 136 ++++++
.../calcMethod/mapNearest/mapNearestMethod.C | 424 ++++++++++++++++++
.../calcMethod/mapNearest/mapNearestMethod.H | 152 +++++++
.../meshToMeshMethod/meshToMeshMethod.C | 274 +++++++++++
.../meshToMeshMethod/meshToMeshMethod.H | 188 ++++++++
.../meshToMeshMethod/meshToMeshMethodI.H | 44 ++
.../meshToMeshMethod/meshToMeshMethodNew.C | 65 +++
.../meshToMeshNew/meshToMeshNew.C | 328 +-------------
.../meshToMeshNew/meshToMeshNew.H | 162 +------
15 files changed, 1920 insertions(+), 969 deletions(-)
delete mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C
delete mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMapNearest.C
rename src/sampling/meshToMeshInterpolation/meshToMeshNew/{calcCellVolumeWeight.C => calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C} (53%)
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.C
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.C
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.C
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
create mode 100644 src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 309cf09d7ab..6d5b7ad19ea 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -60,10 +60,14 @@ $(meshToMesh)/calculateMeshToMeshAddressing.C
$(meshToMesh)/calculateMeshToMeshWeights.C
meshToMeshNew = meshToMeshInterpolation/meshToMeshNew
-$(meshToMeshNew)/calcDirect.C
-$(meshToMeshNew)/calcMapNearest.C
-$(meshToMeshNew)/calcCellVolumeWeight.C
$(meshToMeshNew)/meshToMeshNew.C
$(meshToMeshNew)/meshToMeshNewParallelOps.C
+meshToMeshNewMethods = meshToMeshInterpolation/meshToMeshNew/calcMethod
+$(meshToMeshNewMethods)/meshToMeshMethod/meshToMeshMethod.C
+$(meshToMeshNewMethods)/meshToMeshMethod/meshToMeshMethodNew.C
+$(meshToMeshNewMethods)/cellVolumeWeight/cellVolumeWeightMethod.C
+$(meshToMeshNewMethods)/direct/directMethod.C
+$(meshToMeshNewMethods)/mapNearest/mapNearestMethod.C
+
LIB = $(FOAM_LIBBIN)/libsampling
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C
deleted file mode 100644
index 8e6a17511ef..00000000000
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcDirect.C
+++ /dev/null
@@ -1,159 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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
deleted file mode 100644
index 0276ee56049..00000000000
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMapNearest.C
+++ /dev/null
@@ -1,265 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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/calcCellVolumeWeight.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
similarity index 53%
rename from src/sampling/meshToMeshInterpolation/meshToMeshNew/calcCellVolumeWeight.C
rename to src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
index 9cd8e417f1c..67694afd9d0 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcCellVolumeWeight.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.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) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -23,15 +23,79 @@ License
\*---------------------------------------------------------------------------*/
-#include "meshToMeshNew.H"
-#include "tetOverlapVolume.H"
+#include "cellVolumeWeightMethod.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "addToRunTimeSelectionTable.H"
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-void Foam::meshToMeshNew::calcCellVolumeWeight
+namespace Foam
+{
+ defineTypeNameAndDebug(cellVolumeWeightMethod, 0);
+ addToRunTimeSelectionTable
+ (
+ meshToMeshMethod,
+ cellVolumeWeightMethod,
+ components
+ );
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+bool Foam::cellVolumeWeightMethod::findInitialSeeds
(
- const polyMesh& src,
- const polyMesh& tgt,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+) const
+{
+ const cellList& srcCells = src_.cells();
+ const faceList& srcFaces = src_.faces();
+ const pointField& srcPts = src_.points();
+
+ for (label i = startSeedI; i < srcCellIDs.size(); i++)
+ {
+ label srcI = srcCellIDs[i];
+
+ if (mapFlag[srcI])
+ {
+ const pointField
+ pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
+
+ forAll(pts, ptI)
+ {
+ const point& pt = pts[ptI];
+ label tgtI = tgt_.cellTree().findInside(pt);
+
+ if (tgtI != -1 && intersect(srcI, tgtI))
+ {
+ srcSeedI = srcI;
+ tgtSeedI = tgtI;
+
+ return true;
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "could not find starting seed" << endl;
+ }
+
+ return false;
+}
+
+
+void Foam::cellVolumeWeightMethod::calculateAddressing
+(
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
const label srcSeedI,
const label tgtSeedI,
const labelList& srcCellIDs,
@@ -42,11 +106,11 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
label srcCellI = srcSeedI;
label tgtCellI = tgtSeedI;
- List<DynamicList<label> > srcToTgtAddr(src.nCells());
- List<DynamicList<scalar> > srcToTgtWght(src.nCells());
+ 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<DynamicList<label> > tgtToSrcAddr(tgt_.nCells());
+ List<DynamicList<scalar> > tgtToSrcWght(tgt_.nCells());
// list of tgt cell neighbour cells
DynamicList<label> nbrTgtCells(10);
@@ -55,10 +119,10 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
DynamicList<label> visitedTgtCells(10);
// list to keep track of tgt cells used to seed src cells
- labelList seedCells(src.nCells(), -1);
+ labelList seedCells(src_.nCells(), -1);
seedCells[srcCellI] = tgtCellI;
- const scalarField& srcVol = src.cellVolumes();
+ const scalarField& srcVol = src_.cellVolumes();
do
{
@@ -67,14 +131,14 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
// append initial target cell and neighbours
nbrTgtCells.append(tgtCellI);
- appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+ appendNbrCells(tgtCellI, tgt_, visitedTgtCells, nbrTgtCells);
do
{
tgtCellI = nbrTgtCells.remove();
visitedTgtCells.append(tgtCellI);
- scalar vol = interVol(src, tgt, srcCellI, tgtCellI);
+ scalar vol = interVol(srcCellI, tgtCellI);
// accumulate addressing and weights for valid intersection
if (vol/srcVol[srcCellI] > tolerance_)
@@ -86,7 +150,7 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
tgtToSrcAddr[tgtCellI].append(srcCellI);
tgtToSrcWght[tgtCellI].append(vol);
- appendNbrCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+ appendNbrCells(tgtCellI, tgt_, visitedTgtCells, nbrTgtCells);
// accumulate intersection volume
V_ += vol;
@@ -102,8 +166,6 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
startSeedI,
srcCellI,
tgtCellI,
- src,
- tgt,
srcCellIDs,
mapFlag,
visitedTgtCells,
@@ -113,34 +175,32 @@ void Foam::meshToMeshNew::calcCellVolumeWeight
while (srcCellI != -1);
// transfer addressing into persistent storage
- forAll(srcToTgtCellAddr_, i)
+ forAll(srcToTgtCellAddr, i)
{
- srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]);
- srcToTgtCellWght_[i].transfer(srcToTgtWght[i]);
+ srcToTgtCellAddr[i].transfer(srcToTgtAddr[i]);
+ srcToTgtCellWght[i].transfer(srcToTgtWght[i]);
}
- forAll(tgtToSrcCellAddr_, i)
+ forAll(tgtToSrcCellAddr, i)
{
- tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]);
- tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]);
+ tgtToSrcCellAddr[i].transfer(tgtToSrcAddr[i]);
+ tgtToSrcCellWght[i].transfer(tgtToSrcWght[i]);
}
}
-void Foam::meshToMeshNew::setNextCells
+void Foam::cellVolumeWeightMethod::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];
+ const labelList& srcNbrCells = src_.cellCells()[srcCellI];
// set possible seeds for later use by querying all src cell neighbours
// with all visited target cells
@@ -155,7 +215,7 @@ void Foam::meshToMeshNew::setNextCells
{
label cellT = visitedCells[j];
- if (intersect(src, tgt, cellS, cellT))
+ if (intersect(cellS, cellT))
{
seedCells[cellS] = cellT;
@@ -211,8 +271,6 @@ void Foam::meshToMeshNew::setNextCells
bool restart =
findInitialSeeds
(
- src,
- tgt,
srcCellIDs,
mapFlag,
startSeedI,
@@ -233,68 +291,91 @@ void Foam::meshToMeshNew::setNextCells
}
-bool Foam::meshToMeshNew::intersect
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::cellVolumeWeightMethod::cellVolumeWeightMethod
(
const polyMesh& src,
- const polyMesh& tgt,
- const label srcCellI,
- const label tgtCellI
-) const
-{
- scalar threshold = tolerance_*src.cellVolumes()[srcCellI];
+ const polyMesh& tgt
+)
+:
+ meshToMeshMethod(src, tgt)
+{}
- tetOverlapVolume overlapEngine;
- treeBoundBox bbTgtCell
- (
- pointField
- (
- tgt.points(),
- tgt.cellPoints()[tgtCellI]
- )
- );
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::cellVolumeWeightMethod::~cellVolumeWeightMethod()
+{}
- return overlapEngine.cellCellOverlapMinDecomp
- (
- src,
- srcCellI,
- tgt,
- tgtCellI,
- bbTgtCell,
- threshold
- );
-}
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-Foam::scalar Foam::meshToMeshNew::interVol
+void Foam::cellVolumeWeightMethod::calculate
(
- const polyMesh& src,
- const polyMesh& tgt,
- const label srcCellI,
- const label tgtCellI
-) const
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToSrcAddr,
+ scalarListList& tgtToSrcWght
+)
{
- tetOverlapVolume overlapEngine;
-
- treeBoundBox bbTgtCell
+ bool ok = initialise
(
- pointField
- (
- tgt.points(),
- tgt.cellPoints()[tgtCellI]
- )
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght
);
- scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
- (
- src,
- srcCellI,
- tgt,
- tgtCellI,
- bbTgtCell
- );
+ if (!ok)
+ {
+ return;
+ }
+
+ // (potentially) participating source mesh cells
+ const labelList srcCellIDs(maskCells());
- return vol;
+ // list to keep track of whether src cell can be mapped
+ boolList mapFlag(src_.nCells(), false);
+ UIndirectList<bool>(mapFlag, srcCellIDs) = true;
+
+ // find initial point in tgt mesh
+ label srcSeedI = -1;
+ label tgtSeedI = -1;
+ label startSeedI = 0;
+
+ bool startWalk =
+ findInitialSeeds
+ (
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcSeedI,
+ tgtSeedI
+ );
+
+ if (startWalk)
+ {
+ calculateAddressing
+ (
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght,
+ srcSeedI,
+ tgtSeedI,
+ srcCellIDs,
+ mapFlag,
+ startSeedI
+ );
+ }
+ else
+ {
+ // if meshes are collocated, after inflating the source mesh bounding
+ // box tgt mesh cells may be transferred, but may still not overlap
+ // with the source mesh
+ return;
+ }
}
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
new file mode 100644
index 00000000000..ad38d762465
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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/>.
+
+Class
+ Foam::cellVolumeWeightMethod
+
+Description
+ Cell-volume-weighted mesh-to-mesh interpolation class
+
+ Volume conservative.
+
+SourceFiles
+ cellVolumeWeightMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellVolumeWeightMethod_H
+#define cellVolumeWeightMethod_H
+
+#include "meshToMeshMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class cellVolumeWeightMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellVolumeWeightMethod
+:
+ public meshToMeshMethod
+{
+protected:
+
+ // Protected Member Functions
+
+ //- Find indices of overlapping cells in src and tgt meshes - returns
+ // true if found a matching pair
+ bool findInitialSeeds
+ (
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+ //- Calculate the mesh-to-mesh addressing and weights
+ void calculateAddressing
+ (
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs,
+ boolList& mapFlag,
+ label& startSeedI
+ );
+
+ //- Set the next cells in the advancing front algorithm
+ void setNextCells
+ (
+ label& startSeedI,
+ label& srcCellI,
+ label& tgtCellI,
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const DynamicList<label>& visitedCells,
+ labelList& seedCells
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ cellVolumeWeightMethod(const cellVolumeWeightMethod&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const cellVolumeWeightMethod&);
+
+
+public:
+
+ //- Run-time type information
+ TypeName("cellVolumeWeight");
+
+ //- Construct from source and target meshes
+ cellVolumeWeightMethod(const polyMesh& src, const polyMesh& tgt);
+
+ //- Destructor
+ virtual ~cellVolumeWeightMethod();
+
+
+ // Member Functions
+
+ // Evaluate
+
+ //- Calculate addressing and weights
+ virtual void calculate
+ (
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToTgtAddr,
+ scalarListList& tgtToTgtWght
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.C
new file mode 100644
index 00000000000..9f3674700cd
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.C
@@ -0,0 +1,303 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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 "directMethod.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(directMethod, 0);
+ addToRunTimeSelectionTable(meshToMeshMethod, directMethod, components);
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+bool Foam::directMethod::findInitialSeeds
+(
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+) const
+{
+ const cellList& srcCells = src_.cells();
+ const faceList& srcFaces = src_.faces();
+ const pointField& srcPts = src_.points();
+
+ for (label i = startSeedI; i < srcCellIDs.size(); i++)
+ {
+ label srcI = srcCellIDs[i];
+
+ if (mapFlag[srcI])
+ {
+ const pointField
+ pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
+
+ forAll(pts, ptI)
+ {
+ const point& pt = pts[ptI];
+ label tgtI = tgt_.cellTree().findInside(pt);
+
+ if (tgtI != -1 && intersect(srcI, tgtI))
+ {
+ srcSeedI = srcI;
+ tgtSeedI = tgtI;
+
+ return true;
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "could not find starting seed" << endl;
+ }
+
+ return false;
+}
+
+
+void Foam::directMethod::calculateAddressing
+(
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs, // not used
+ boolList& mapFlag,
+ label& startSeedI
+)
+{
+ // store a list of src cells already mapped
+ labelList srcTgtSeed(src_.nCells(), -1);
+
+ List<DynamicList<label> > srcToTgt(src_.nCells());
+ List<DynamicList<label> > tgtToSrc(tgt_.nCells());
+
+ DynamicList<label> srcSeeds(10);
+
+ 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
+ mapFlag[srcCellI] = false;
+
+ // accumulate intersection volume
+ V_ += srcVc[srcCellI];
+
+ // find new source seed cell
+ appendToDirectSeeds
+ (
+ mapFlag,
+ 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(1, v);
+ }
+
+ forAll(tgtToSrcCellAddr, i)
+ {
+ scalar v = tgtVc[i];
+ tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
+ tgtToSrcCellWght[i] = scalarList(1, v);
+ }
+}
+
+
+void Foam::directMethod::appendToDirectSeeds
+(
+ 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;
+ }
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::directMethod::directMethod
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+)
+:
+ meshToMeshMethod(src, tgt)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::directMethod::~directMethod()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::directMethod::calculate
+(
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToSrcAddr,
+ scalarListList& tgtToSrcWght
+)
+{
+ bool ok = initialise
+ (
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght
+ );
+
+ if (!ok)
+ {
+ return;
+ }
+
+ // (potentially) participating source mesh cells
+ const labelList srcCellIDs(maskCells());
+
+ // list to keep track of whether src cell can be mapped
+ boolList mapFlag(src_.nCells(), false);
+ UIndirectList<bool>(mapFlag, srcCellIDs) = true;
+
+ // find initial point in tgt mesh
+ label srcSeedI = -1;
+ label tgtSeedI = -1;
+ label startSeedI = 0;
+
+ bool startWalk =
+ findInitialSeeds
+ (
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcSeedI,
+ tgtSeedI
+ );
+
+ if (startWalk)
+ {
+ calculateAddressing
+ (
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght,
+ srcSeedI,
+ tgtSeedI,
+ srcCellIDs,
+ mapFlag,
+ startSeedI
+ );
+ }
+ else
+ {
+ // if meshes are collocated, after inflating the source mesh bounding
+ // box tgt mesh cells may be transferred, but may still not overlap
+ // with the source mesh
+ return;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.H
new file mode 100644
index 00000000000..1c9f489a83c
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/direct/directMethod.H
@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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/>.
+
+Class
+ Foam::directMethod
+
+Description
+ Direct (one-to-one cell correspondence) mesh-to-mesh interpolation class
+
+ Volume conservative.
+
+SourceFiles
+ directMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMethod_H
+#define directMethod_H
+
+#include "meshToMeshMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class directMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMethod
+:
+ public meshToMeshMethod
+{
+protected:
+
+ // Protected Member Functions
+
+ //- Find indices of overlapping cells in src and tgt meshes - returns
+ // true if found a matching pair
+ bool findInitialSeeds
+ (
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+ //- Calculate the mesh-to-mesh addressing and weights
+ void calculateAddressing
+ (
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs,
+ boolList& mapFlag,
+ label& startSeedI
+ );
+
+ //- Append to list of src mesh seed indices
+ void appendToDirectSeeds
+ (
+ boolList& mapFlag,
+ labelList& srcTgtSeed,
+ DynamicList<label>& srcSeeds,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ directMethod(const directMethod&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const directMethod&);
+
+
+public:
+
+ //- Run-time type information
+ TypeName("direct");
+
+ //- Construct from source and target meshes
+ directMethod(const polyMesh& src, const polyMesh& tgt);
+
+ //- Destructor
+ virtual ~directMethod();
+
+
+ // Member Functions
+
+ // Evaluate
+
+ //- Calculate addressing and weights
+ virtual void calculate
+ (
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToTgtAddr,
+ scalarListList& tgtToTgtWght
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.C
new file mode 100644
index 00000000000..9327c3abb19
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.C
@@ -0,0 +1,424 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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 "mapNearestMethod.H"
+#include "pointIndexHit.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(mapNearestMethod, 0);
+ addToRunTimeSelectionTable(meshToMeshMethod, mapNearestMethod, components);
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+bool Foam::mapNearestMethod::findInitialSeeds
+(
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+) const
+{
+ const cellList& srcCells = src_.cells();
+ const faceList& srcFaces = src_.faces();
+ const pointField& srcPts = src_.points();
+
+ for (label i = startSeedI; i < srcCellIDs.size(); i++)
+ {
+ label srcI = srcCellIDs[i];
+
+ if (mapFlag[srcI])
+ {
+ const pointField
+ pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
+
+ const point& pt = pts[0];
+ pointIndexHit hit = tgt_.cellTree().findNearest(pt, GREAT);
+
+ if (hit.hit())
+ {
+ srcSeedI = srcI;
+ tgtSeedI = hit.index();
+
+ return true;
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "bool Foam::mapNearestMethod::findInitialSeeds"
+ "("
+ "const labelList&, "
+ "const boolList&, "
+ "const label, "
+ "label&, "
+ "label&"
+ ") const"
+ )
+ << "Unable to find nearest target cell"
+ << " for source cell " << srcI
+ << " with centre "
+ << srcCells[srcI].centre(srcPts, srcFaces)
+ << abort(FatalError);
+ }
+
+ break;
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "could not find starting seed" << endl;
+ }
+
+ return false;
+}
+
+
+void Foam::mapNearestMethod::calculateAddressing
+(
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
+ 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,
+ 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(tgtCellI, tgtToSrc);
+
+ findNearestCell(tgt_, src_, tgtCellI, srcCellI);
+
+ tgtToSrc[tgtCellI].append(srcCellI);
+ }
+ }
+
+ // transfer addressing into persistent storage
+ forAll(srcToTgtCellAddr, i)
+ {
+ scalar v = srcVc[i];
+ srcToTgtCellAddr[i].transfer(srcToTgt[i]);
+ srcToTgtCellWght[i] = scalarList(1, v);
+ }
+
+ forAll(tgtToSrcCellAddr, i)
+ {
+ scalar v = tgtVc[i];
+ tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
+ tgtToSrcCellWght[i] = scalarList(1, v);
+ }
+}
+
+
+void Foam::mapNearestMethod::findNearestCell
+(
+ const polyMesh& mesh1,
+ const polyMesh& mesh2,
+ const label cell1,
+ label& cell2
+) const
+{
+ const vectorField& Cc1 = mesh1.cellCentres();
+ const vectorField& Cc2 = mesh2.cellCentres();
+
+ const vector& p1 = Cc1[cell1];
+
+ DynamicList<label> cells2(10);
+ cells2.append(cell2);
+
+ DynamicList<label> visitedCells(10);
+
+ scalar d = GREAT;
+
+ do
+ {
+ label c2 = cells2.remove();
+ visitedCells.append(c2);
+
+ scalar dTest = magSqr(Cc2[c2] - p1);
+ if (dTest < d)
+ {
+ cell2 = c2;
+ d = dTest;
+ appendNbrCells(cell2, mesh2, visitedCells, cells2);
+ }
+
+ } while (cells2.size() > 0);
+}
+
+
+void Foam::mapNearestMethod::setNextNearestCells
+(
+ label& startSeedI,
+ label& srcCellI,
+ label& tgtCellI,
+ boolList& mapFlag,
+ const labelList& srcCellIDs
+) const
+{
+ 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
+ (
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcCellI,
+ tgtCellI
+ );
+}
+
+
+Foam::label Foam::mapNearestMethod::findMappedSrcCell
+(
+ 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;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::mapNearestMethod::mapNearestMethod
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+)
+:
+ meshToMeshMethod(src, tgt)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::mapNearestMethod::~mapNearestMethod()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::mapNearestMethod::calculate
+(
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToSrcAddr,
+ scalarListList& tgtToSrcWght
+)
+{
+ bool ok = initialise
+ (
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght
+ );
+
+ if (!ok)
+ {
+ return;
+ }
+
+ // (potentially) participating source mesh cells
+ const labelList srcCellIDs(maskCells());
+
+ // list to keep track of whether src cell can be mapped
+ boolList mapFlag(src_.nCells(), false);
+ UIndirectList<bool>(mapFlag, srcCellIDs) = true;
+
+ // find initial point in tgt mesh
+ label srcSeedI = -1;
+ label tgtSeedI = -1;
+ label startSeedI = 0;
+
+ bool startWalk =
+ findInitialSeeds
+ (
+ srcCellIDs,
+ mapFlag,
+ startSeedI,
+ srcSeedI,
+ tgtSeedI
+ );
+
+ if (startWalk)
+ {
+ calculateAddressing
+ (
+ srcToTgtAddr,
+ srcToTgtWght,
+ tgtToSrcAddr,
+ tgtToSrcWght,
+ srcSeedI,
+ tgtSeedI,
+ srcCellIDs,
+ mapFlag,
+ startSeedI
+ );
+ }
+ else
+ {
+ // if meshes are collocated, after inflating the source mesh bounding
+ // box tgt mesh cells may be transferred, but may still not overlap
+ // with the source mesh
+ return;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.H
new file mode 100644
index 00000000000..813c46683ff
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/mapNearest/mapNearestMethod.H
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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/>.
+
+Class
+ Foam::mapNearestMethod
+
+Description
+ Map nearest mesh-to-mesh interpolation class
+
+ Not volume conservative.
+
+SourceFiles
+ mapNearestMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mapNearestMethod_H
+#define mapNearestMethod_H
+
+#include "meshToMeshMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class mapNearestMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class mapNearestMethod
+:
+ public meshToMeshMethod
+{
+protected:
+
+ // Protected Member Functions
+
+ //- Find indices of overlapping cells in src and tgt meshes - returns
+ // true if found a matching pair
+ bool findInitialSeeds
+ (
+ const labelList& srcCellIDs,
+ const boolList& mapFlag,
+ const label startSeedI,
+ label& srcSeedI,
+ label& tgtSeedI
+ ) const;
+
+ //- Calculate the mesh-to-mesh addressing and weights
+ void calculateAddressing
+ (
+ labelListList& srcToTgtCellAddr,
+ scalarListList& srcToTgtCellWght,
+ labelListList& tgtToSrcCellAddr,
+ scalarListList& tgtToSrcCellWght,
+ const label srcSeedI,
+ const label tgtSeedI,
+ const labelList& srcCellIDs,
+ boolList& mapFlag,
+ label& startSeedI
+ );
+
+ //- Find the nearest cell on mesh2 for cell1 on mesh1
+ void findNearestCell
+ (
+ const polyMesh& mesh1,
+ const polyMesh& mesh2,
+ const label cell1,
+ label& cell2
+ ) const;
+
+ //- Set the next cells for the marching front algorithm
+ void setNextNearestCells
+ (
+ label& startSeedI,
+ label& srcCellI,
+ label& tgtCellI,
+ boolList& mapFlag,
+ const labelList& srcCellIDs
+ ) const;
+
+ //- Find a source cell mapped to target cell tgtCellI
+ label findMappedSrcCell
+ (
+ const label tgtCellI,
+ const List<DynamicList<label> >& tgtToSrc
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ mapNearestMethod(const mapNearestMethod&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const mapNearestMethod&);
+
+
+public:
+
+ //- Run-time type information
+ TypeName("mapNearest");
+
+ //- Construct from source and target meshes
+ mapNearestMethod(const polyMesh& src, const polyMesh& tgt);
+
+ //- Destructor
+ virtual ~mapNearestMethod();
+
+
+ // Member Functions
+
+ // Evaluate
+
+ //- Calculate addressing and weights
+ virtual void calculate
+ (
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToTgtAddr,
+ scalarListList& tgtToTgtWght
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.C
new file mode 100644
index 00000000000..29e066dd542
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.C
@@ -0,0 +1,274 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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 "meshToMeshMethod.H"
+#include "tetOverlapVolume.H"
+#include "OFstream.H"
+#include "Time.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(meshToMeshMethod, 0);
+ defineRunTimeSelectionTable(meshToMeshMethod, components);
+}
+
+Foam::scalar Foam::meshToMeshMethod::tolerance_ = 1e-6;
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::labelList Foam::meshToMeshMethod::maskCells() const
+{
+ boundBox intersectBb
+ (
+ max(src_.bounds().min(), tgt_.bounds().min()),
+ min(src_.bounds().max(), tgt_.bounds().max())
+ );
+
+ intersectBb.inflate(0.01);
+
+ const cellList& srcCells = src_.cells();
+ const faceList& srcFaces = src_.faces();
+ const pointField& srcPts = src_.points();
+
+ DynamicList<label> cells(src_.nCells());
+ forAll(srcCells, srcI)
+ {
+ boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
+ if (intersectBb.overlaps(cellBb))
+ {
+ cells.append(srcI);
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "participating source mesh cells: " << cells.size() << endl;
+ }
+
+ return cells;
+}
+
+
+bool Foam::meshToMeshMethod::intersect
+(
+ 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::meshToMeshMethod::interVol
+(
+ 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::meshToMeshMethod::appendNbrCells
+(
+ const label cellI,
+ const polyMesh& mesh,
+ const DynamicList<label>& visitedCells,
+ DynamicList<label>& nbrCellIDs
+) const
+{
+ const labelList& nbrCells = mesh.cellCells()[cellI];
+
+ // filter out cells already visited from cell neighbours
+ forAll(nbrCells, i)
+ {
+ label nbrCellI = nbrCells[i];
+
+ if
+ (
+ (findIndex(visitedCells, nbrCellI) == -1)
+ && (findIndex(nbrCellIDs, nbrCellI) == -1)
+ )
+ {
+ nbrCellIDs.append(nbrCellI);
+ }
+ }
+}
+
+
+bool Foam::meshToMeshMethod::initialise
+(
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToSrcAddr,
+ scalarListList& tgtToSrcWght
+) const
+{
+ srcToTgtAddr.setSize(src_.nCells());
+ srcToTgtWght.setSize(src_.nCells());
+ tgtToSrcAddr.setSize(tgt_.nCells());
+ tgtToSrcWght.setSize(tgt_.nCells());
+
+ if (!src_.nCells())
+ {
+ return false;
+ }
+ else if (!tgt_.nCells())
+ {
+ if (debug)
+ {
+ Pout<< "mesh interpolation: hhave " << src_.nCells() << " source "
+ << " cells but no target cells" << endl;
+ }
+
+ return false;
+ }
+
+ return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::meshToMeshMethod::meshToMeshMethod
+(
+ const polyMesh& src,
+ const polyMesh& tgt
+)
+:
+ src_(src),
+ tgt_(tgt),
+ V_(0.0)
+{
+ if (!src_.nCells() || !tgt_.nCells())
+ {
+ if (debug)
+ {
+ Pout<< "mesh interpolation: cells not on processor: Source cells = "
+ << src_.nCells() << ", target cells = " << tgt_.nCells()
+ << endl;
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::meshToMeshMethod::~meshToMeshMethod()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::meshToMeshMethod::writeConnectivity
+(
+ const polyMesh& mesh1,
+ const polyMesh& mesh2,
+ const labelListList& mesh1ToMesh2Addr
+) const
+{
+ Pout<< "Source size = " << mesh1.nCells() << endl;
+ Pout<< "Target size = " << mesh2.nCells() << endl;
+
+ word fName("addressing_" + mesh1.name() + "_to_" + mesh2.name());
+
+ if (Pstream::parRun())
+ {
+ fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
+ }
+
+ OFstream os(src_.time().path()/fName + ".obj");
+
+ label vertI = 0;
+ forAll(mesh1ToMesh2Addr, i)
+ {
+ const labelList& addr = mesh1ToMesh2Addr[i];
+ forAll(addr, j)
+ {
+ label cellI = addr[j];
+ const vector& c0 = mesh1.cellCentres()[i];
+
+ const cell& c = mesh2.cells()[cellI];
+ const pointField pts(c.points(mesh2.faces(), mesh2.points()));
+ forAll(pts, j)
+ {
+ const point& p = pts[j];
+ os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+ vertI++;
+ os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
+ << nl;
+ vertI++;
+ os << "l " << vertI - 1 << ' ' << vertI << nl;
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.H
new file mode 100644
index 00000000000..3338fa53495
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethod.H
@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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/>.
+
+Class
+ Foam::meshToMeshMethod
+
+Description
+ Base class for mesh-to-mesh calculation methods
+
+SourceFiles
+ meshToMeshMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshToMeshMethod_H
+#define meshToMeshMethod_H
+
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class meshToMeshMethod Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshToMeshMethod
+{
+
+protected:
+
+ // Protected data
+
+ //- Reference to the source mesh
+ const polyMesh& src_;
+
+ //- Reference to the target mesh
+ const polyMesh& tgt_;
+
+ //- Cell total volume in overlap region [m3]
+ scalar V_;
+
+ //- Tolerance used in volume overlap calculations
+ static scalar tolerance_;
+
+
+ // Protected Member Functions
+
+ //- Return src cell IDs for the overlap region
+ labelList maskCells() const;
+
+ //- Return the true if cells intersect
+ bool intersect(const label srcCellI, const label tgtCellI) const;
+
+ //- Return the intersection volume between two cells
+ scalar interVol(const label srcCellI, const label tgtCellI) const;
+
+ //- Append target cell neihgbour cells to cellIDs list
+ void appendNbrCells
+ (
+ const label tgtCellI,
+ const polyMesh& mesh,
+ const DynamicList<label>& visitedTgtCells,
+ DynamicList<label>& nbrTgtCellIDs
+ ) const;
+
+ bool initialise
+ (
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToTgtAddr,
+ scalarListList& tgtToTgtWght
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ meshToMeshMethod(const meshToMeshMethod&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const meshToMeshMethod&);
+
+
+public:
+
+ //- Run-time type information
+ TypeName("meshToMeshMethod");
+
+ //- Declare runtime constructor selection table
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ meshToMeshMethod,
+ components,
+ (
+ const polyMesh& src,
+ const polyMesh& tgt
+ ),
+ (src, tgt)
+ );
+
+ //- Construct from source and target meshes
+ meshToMeshMethod(const polyMesh& src, const polyMesh& tgt);
+
+ //- Selector
+ static autoPtr<meshToMeshMethod> New
+ (
+ const word& methodName,
+ const polyMesh& src,
+ const polyMesh& tgt
+ );
+
+
+ //- Destructor
+ virtual ~meshToMeshMethod();
+
+
+ // Member Functions
+
+ // Evaluate
+
+ //- Calculate addressing and weights
+ virtual void calculate
+ (
+ labelListList& srcToTgtAddr,
+ scalarListList& srcToTgtWght,
+ labelListList& tgtToTgtAddr,
+ scalarListList& tgtToTgtWght
+ ) = 0;
+
+
+ // Access
+
+ //- Return const access to the source mesh
+ inline const polyMesh& src() const;
+
+ //- Return const access to the target mesh
+ inline const polyMesh& tgt() const;
+
+ //- Return const access to the overlap volume
+ inline scalar V() const;
+
+
+ // Check
+
+ //- Write the connectivity (debugging)
+ void writeConnectivity
+ (
+ const polyMesh& mesh1,
+ const polyMesh& mesh2,
+ const labelListList& mesh1ToMesh2Addr
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "meshToMeshMethodI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodI.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
new file mode 100644
index 00000000000..a1f98a9123c
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+const Foam::polyMesh& Foam::meshToMeshMethod::src() const
+{
+ return src_;
+}
+
+
+const Foam::polyMesh& Foam::meshToMeshMethod::tgt() const
+{
+ return tgt_;
+}
+
+
+Foam::scalar Foam::meshToMeshMethod::V() const
+{
+ return V_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
new file mode 100644
index 00000000000..befb3e6debd
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 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 "meshToMeshMethod.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::meshToMeshMethod> Foam::meshToMeshMethod::New
+(
+ const word& methodName,
+ const polyMesh& src,
+ const polyMesh& tgt
+)
+{
+ if (debug)
+ {
+ Info<< "Selecting AMIMethod " << methodName << endl;
+ }
+
+ componentsConstructorTable::iterator cstrIter =
+ componentsConstructorTablePtr_->find(methodName);
+
+ if (cstrIter == componentsConstructorTablePtr_->end())
+ {
+ FatalErrorIn
+ (
+ "Foam::autoPtr<Foam::meshToMeshMethod> Foam::meshToMeshMethod::New"
+ "("
+ "const word&, "
+ "const polyMesh&, "
+ "const polyMesh&"
+ ")"
+ ) << "Unknown meshToMesh type "
+ << methodName << nl << nl
+ << "Valid meshToMesh types are:" << nl
+ << componentsConstructorTablePtr_->sortedToc() << exit(FatalError);
+ }
+
+ return autoPtr<meshToMeshMethod>(cstrIter()(src, tgt));
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
index 647b78b6cd6..27f0bd592c3 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
@@ -24,14 +24,9 @@ License
\*---------------------------------------------------------------------------*/
#include "meshToMeshNew.H"
-#include "OFstream.H"
#include "Time.H"
#include "globalIndex.H"
-#include "mergePoints.H"
-#include "treeBoundBox.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
-#include "ListOps.H"
+#include "meshToMeshMethod.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -55,55 +50,9 @@ namespace Foam
meshToMeshNew::interpolationMethodNames_;
}
-Foam::scalar Foam::meshToMeshNew::tolerance_ = 1e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-void Foam::meshToMeshNew::writeConnectivity
-(
- const polyMesh& src,
- const polyMesh& tgt,
- const labelListList& srcToTargetAddr
-) const
-{
- Pout<< "Source size = " << src.nCells() << endl;
- Pout<< "Target size = " << tgt.nCells() << endl;
-
- word fName("addressing_" + src.name() + "_to_" + tgt.name());
-
- if (Pstream::parRun())
- {
- fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
- }
-
- OFstream os(src.time().path()/fName + ".obj");
-
- label vertI = 0;
- forAll(srcToTargetAddr, i)
- {
- const labelList& tgtAddress = srcToTargetAddr[i];
- forAll(tgtAddress, j)
- {
- label tgtI = tgtAddress[j];
- const vector& c0 = src.cellCentres()[i];
-
- const cell& c = tgt.cells()[tgtI];
- const pointField pts(c.points(tgt.faces(), tgt.points()));
- forAll(pts, j)
- {
- const point& p = pts[j];
- os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
- vertI++;
- os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
- << nl;
- vertI++;
- os << "l " << vertI - 1 << ' ' << vertI << nl;
- }
- }
- }
-}
-
-
Foam::labelList Foam::meshToMeshNew::maskCells
(
const polyMesh& src,
@@ -141,146 +90,6 @@ Foam::labelList Foam::meshToMeshNew::maskCells
}
-bool Foam::meshToMeshNew::findInitialSeeds
-(
- const polyMesh& src,
- const polyMesh& tgt,
- const labelList& srcCellIDs,
- const boolList& mapFlag,
- const label startSeedI,
- label& srcSeedI,
- label& tgtSeedI
-) const
-{
- const cellList& srcCells = src.cells();
- const faceList& srcFaces = src.faces();
- const pointField& srcPts = src.points();
-
- for (label i = startSeedI; i < srcCellIDs.size(); i++)
- {
- label srcI = srcCellIDs[i];
-
- if (mapFlag[srcI])
- {
- const pointField
- pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
-
- switch (method_)
- {
- case imDirect:
- case imCellVolumeWeight:
- {
- forAll(pts, ptI)
- {
- const point& pt = pts[ptI];
- label tgtI = tgt.cellTree().findInside(pt);
-
- if (tgtI != -1 && intersect(src, tgt, srcI, tgtI))
- {
- srcSeedI = srcI;
- tgtSeedI = tgtI;
-
- return true;
- }
- }
-
- break;
- }
- case imMapNearest:
- {
- const point& pt = pts[0];
- pointIndexHit hit = tgt.cellTree().findNearest(pt, GREAT);
-
- if (hit.hit())
- {
- srcSeedI = srcI;
- tgtSeedI = hit.index();
-
- return true;
- }
- else
- {
- FatalErrorIn
- (
- "bool Foam::meshToMeshNew::findInitialSeeds"
- "("
- "const polyMesh&, "
- "const polyMesh&, "
- "const labelList&, "
- "const boolList&, "
- "const label, "
- "label&, "
- "label&"
- ") const"
- )
- << "Unable to find nearest target cell"
- << " for source cell " << srcI
- << " with centre "
- << srcCells[srcI].centre(srcPts, srcFaces)
- << abort(FatalError);
- }
-
- break;
- }
- default:
- {
- FatalErrorIn
- (
- "bool Foam::meshToMeshNew::findInitialSeeds"
- "("
- "const polyMesh&, "
- "const polyMesh&, "
- "const labelList&, "
- "const boolList&, "
- "const label, "
- "label&, "
- "label&"
- ") const"
- )
- << "Unhandled method: "
- << interpolationMethodNames_[method_]
- << abort(FatalError);
- }
- }
- }
- }
-
- if (debug)
- {
- Pout<< "could not find starting seed" << endl;
- }
-
- return false;
-}
-
-
-void Foam::meshToMeshNew::appendNbrCells
-(
- const label cellI,
- const polyMesh& mesh,
- const DynamicList<label>& visitedCells,
- DynamicList<label>& nbrCellIDs
-) const
-{
- const labelList& nbrCells = mesh.cellCells()[cellI];
-
- // filter out cells already visited from cell neighbours
- forAll(nbrCells, i)
- {
- label nbrCellI = nbrCells[i];
-
- if
- (
- (findIndex(visitedCells, nbrCellI) == -1)
- && (findIndex(nbrCellIDs, nbrCellI) == -1)
- )
- {
- nbrCellIDs.append(nbrCellI);
- }
- }
-}
-
-
void Foam::meshToMeshNew::normaliseWeights
(
const word& descriptor,
@@ -324,124 +133,29 @@ void Foam::meshToMeshNew::calcAddressing
const polyMesh& tgt
)
{
- srcToTgtCellAddr_.setSize(src.nCells());
- srcToTgtCellWght_.setSize(src.nCells());
-
- tgtToSrcCellAddr_.setSize(tgt.nCells());
- tgtToSrcCellWght_.setSize(tgt.nCells());
-
- if (!src.nCells() || !tgt.nCells())
- {
- if (debug)
- {
- Pout<< "mesh interpolation: cells not on processor: Source cells = "
- << src.nCells() << ", target cells = " << tgt.nCells()
- << endl;
- }
- }
-
- if (!src.nCells())
- {
- return;
- }
- else if (!tgt.nCells())
- {
- if (debug)
- {
- Pout<< "mesh interpolation: hhave " << src.nCells() << " source "
- << " cells but no target cells" << endl;
- }
-
- return;
- }
-
- // (potentially) participating source mesh cells
- const labelList srcCellIDs = maskCells(src, tgt);
-
- // list to keep track of whether src cell can be mapped
- boolList mapFlag(src.nCells(), false);
- UIndirectList<bool>(mapFlag, srcCellIDs) = true;
-
- // find initial point in tgt mesh
- label srcSeedI = -1;
- label tgtSeedI = -1;
- label startSeedI = 0;
-
- bool startWalk =
- findInitialSeeds
+ autoPtr<meshToMeshMethod> methodPtr
+ (
+ meshToMeshMethod::New
(
+ interpolationMethodNames_[method_],
src,
- tgt,
- srcCellIDs,
- mapFlag,
- startSeedI,
- srcSeedI,
- tgtSeedI
- );
-
- if (!startWalk)
- {
- // if meshes are collocated, after inflating the source mesh bounding
- // box tgt mesh cells may be transferred, but may still not overlap
- // with the source mesh
- return;
- }
-
+ tgt
+ )
+ );
- switch (method_)
- {
- case imDirect:
- {
- calcDirect(src, tgt, srcSeedI, tgtSeedI);
- break;
- }
- case imMapNearest:
- {
- calcMapNearest
- (
- src,
- tgt,
- srcSeedI,
- tgtSeedI,
- srcCellIDs,
- mapFlag,
- startSeedI
- );
- break;
- }
- case imCellVolumeWeight:
- {
- calcCellVolumeWeight
- (
- src,
- tgt,
- srcSeedI,
- tgtSeedI,
- srcCellIDs,
- mapFlag,
- startSeedI
- );
- break;
- }
- default:
- {
- FatalErrorIn
- (
- "void Foam::meshToMeshNew::calcAddressing"
- "("
- "const polyMesh&, "
- "const polyMesh&"
- ")"
- )
- << "Unknown interpolation method"
- << abort(FatalError);
- }
- }
+ methodPtr->calculate
+ (
+ srcToTgtCellAddr_,
+ srcToTgtCellWght_,
+ tgtToSrcCellAddr_,
+ tgtToSrcCellWght_
+ );
+ V_ = methodPtr->V();
if (debug)
{
- writeConnectivity(src, tgt, srcToTgtCellAddr_);
+ methodPtr->writeConnectivity(src, tgt, srcToTgtCellAddr_);
}
}
@@ -688,6 +402,12 @@ Foam::meshToMeshNew::patchAMIs() const
{
if (patchAMIs_.empty())
{
+ const word amiMethod =
+ AMIPatchToPatchInterpolation::interpolationMethodToWord
+ (
+ interpolationMethodAMI(method_)
+ );
+
patchAMIs_.setSize(srcPatchID_.size());
forAll(srcPatchID_, i)
@@ -700,7 +420,7 @@ Foam::meshToMeshNew::patchAMIs() const
Info<< "Creating AMI between source patch " << srcPP.name()
<< " and target patch " << tgtPP.name()
- << " using " << interpolationMethodAMI(method_)
+ << " using " << amiMethod
<< endl;
Info<< incrIndent;
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
index 19ec3b3e0d6..ab05a058ae9 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
@@ -27,20 +27,14 @@ 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.
+ Mapping is performed using a run-time selectable interpolation mothod
+SeeAlso
+ meshToMeshMethod
SourceFiles
- calcDirect.C
- calcMapNearest.C
- calcCellVolumeWeight.C
meshToMeshNew.C
+ meshToMeshNewParallelOps.C
meshToMeshNewTemplates.C
\*---------------------------------------------------------------------------*/
@@ -128,9 +122,6 @@ private:
//- Target map pointer - parallel running only
autoPtr<mapDistribute> tgtMapPtr_;
- //- Tolerance used in volume overlap calculations
- static scalar tolerance_;
-
// Private Member Functions
@@ -138,39 +129,9 @@ private:
template<class Type>
void add(UList<Type>& fld, const label offset) const;
- //- Write the connectivity (debugging)
- void writeConnectivity
- (
- const polyMesh& src,
- const polyMesh& tgt,
- const labelListList& srcToTargetAddr
- ) const;
-
//- Return src cell IDs for the overlap region
labelList maskCells(const polyMesh& src, const polyMesh& tgt) const;
- //- Find indices of overlapping cells in src and tgt meshes - returns
- // true if found a matching pair
- bool findInitialSeeds
- (
- const polyMesh& src,
- const polyMesh& tgt,
- const labelList& srcCellIDs,
- const boolList& mapFlag,
- const label startSeedI,
- label& srcSeedI,
- 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
(
@@ -180,121 +141,6 @@ private:
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
- (
- const polyMesh& src,
- const polyMesh& tgt,
- boolList& mapFlag,
- labelList& srcTgtSeed,
- DynamicList<label>& srcSeeds,
- label& srcSeedI,
- label& tgtSeedI
- ) const;
-
- // 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 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
- );
-
- //- 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
- );
-
- //- Find source cell for target cell
- label findMappedSrcCell
- (
- const polyMesh& tgt,
- const label tgtCellI,
- const List<DynamicList<label> >& tgtToSrc
- ) const;
-
-
- // Cell volume weighted (non-conformal) interpolation
-
- //- Main driver routine for cell volume weighted interpolation
- void calcCellVolumeWeight
- (
- const polyMesh& src,
- const polyMesh& tgt,
- const label srcSeedI,
- const label tgtSeedI,
- const labelList& srcCellIDs,
- boolList& mapFlag,
- label& startSeedI
- );
-
- //- Set the next cells in the advancing front algorithm
- void 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;
-
- //- Return the true if cells intersect
- bool intersect
- (
- const polyMesh& src,
- const polyMesh& tgt,
- const label srcCellI,
- const label tgtCellI
- ) const;
-
- //- Return the intersection volume between two cells
- scalar interVol
- (
- const polyMesh& src,
- const polyMesh& tgt,
- const label srcCellI,
- const label tgtCellI
- ) const;
-
-
//- Calculate the addressing between overalping regions of src and tgt
// meshes
void calcAddressing(const polyMesh& src, const polyMesh& tgt);
--
GitLab