Commit 6000d2fb authored by andy's avatar andy
Browse files

ENH: Refectored mesh-to-mesh interpolation methods - now run-time selectable

parent 587ca3ee
......@@ -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
......@@ -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;
}
}
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
// ************************************************************************* //
......@@ -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,29 +23,91 @@ License
\*---------------------------------------------------------------------------*/
#include "meshToMeshNew.H"
#include "directMethod.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
void Foam::meshToMeshNew::calcDirect
namespace Foam
{
defineTypeNameAndDebug(directMethod, 0);
addToRunTimeSelectionTable(meshToMeshMethod, directMethod, components);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
bool Foam::directMethod::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)
{