Commit f4d592c3 authored by andy's avatar andy
Browse files

ENH: Added map nearest option to mesh-to-mesh interpolation

parent 5baa8f5a
......@@ -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
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
boolList tgtMapFlag(tgt.nCells(), true);
do
{
srcToTgt[srcCellI].append(tgtCellI);
// find nearest tgt cell
findNearestCell(src, tgt, srcCellI, tgtCellI);
// store src/tgt cell pair
tgtToSrc[tgtCellI].append(srcCellI);
// mark source cell srcCellI and tgtCellI as matched
mapFlag[srcCellI] = false;
tgtMapFlag[tgtCellI] = false;
// accumulate intersection volume
V_ += srcVc[srcCellI];
// find new source cell
setNextNearestCells
(
startSeedI,
srcCellI,
tgtCellI,
mapFlag,
src,
tgt,
srcCellIDs
);
}
while (srcCellI >= 0);
// If there are more target cells than source cells, some target cells
// will not yet be mapped
forAll(tgtMapFlag, tgtCellI)
{
if (tgtMapFlag[tgtCellI])
{
label srcCellI = findMappedSrcCell(tgt, tgtCellI, tgtToSrc);
findNearestCell(tgt, src, tgtCellI, srcCellI);
tgtToSrc[tgtCellI].append(srcCellI);
}
}
// transfer addressing into persistent storage
// note: always 1 target cell per source cell (srcToTgt)
// can be multiple source cells per target cell (tgtToSrc)
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];
scalarList w(tgtToSrc[i].size(), v);
forAll(w, j)
{
w[j] /= w.size();
}
tgtToSrcCellWght_[i] = scalarList(w, v);
tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]);
}
}
void Foam::meshToMeshNew::findNearestCell
(
const polyMesh& src,
const polyMesh& tgt,
const label srcCellI,
label& tgtCellI
)