ENH: Updated meshToMesh to be able to use the new processorLOD

This commit is contained in:
Andrew Heather
2018-06-20 07:57:25 +01:00
parent d5422239f3
commit 1c64911367
7 changed files with 281 additions and 146 deletions

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -46,6 +46,7 @@ void mapConsistentMesh
const fvMesh& meshTarget,
const word& mapMethod,
const word& AMIMapMethod,
const word& procMapMethod,
const bool subtract,
const wordHashSet& selectedFields,
const bool noLagrangian
@ -54,7 +55,14 @@ void mapConsistentMesh
Info<< nl << "Consistently creating and mapping fields for time "
<< meshSource.time().timeName() << nl << endl;
meshToMesh interp(meshSource, meshTarget, mapMethod, AMIMapMethod);
meshToMesh interp
(
meshSource,
meshTarget,
mapMethod,
AMIMapMethod,
meshToMesh::procMapMethodNames_[procMapMethod]
);
if (subtract)
{
@ -85,6 +93,7 @@ void mapSubMesh
const wordList& cuttingPatches,
const word& mapMethod,
const word& AMIMapMethod,
const word& procMapMethod,
const bool subtract,
const wordHashSet& selectedFields,
const bool noLagrangian
@ -100,7 +109,8 @@ void mapSubMesh
mapMethod,
AMIMapMethod,
patchMap,
cuttingPatches
cuttingPatches,
meshToMesh::procMapMethodNames_[procMapMethod]
);
if (subtract)
@ -171,6 +181,12 @@ int main(int argc, char *argv[])
"word",
"specify the patch mapping method (direct|mapNearest|faceAreaWeight)"
);
argList::addOption
(
"procMapMethod",
"word",
"specify the processor distribution map method (AABB|LOD)"
);
argList::addBoolOption
(
"subtract",
@ -217,7 +233,7 @@ int main(int argc, char *argv[])
word mapMethod = meshToMesh::interpolationMethodNames_
[
meshToMesh::imCellVolumeWeight
meshToMesh::interpolationMethod::imCellVolumeWeight
];
if (args.readIfPresent("mapMethod", mapMethod))
@ -225,7 +241,6 @@ int main(int argc, char *argv[])
Info<< "Mapping method: " << mapMethod << endl;
}
word patchMapMethod;
if (meshToMesh::interpolationMethodNames_.found(mapMethod))
{
@ -233,12 +248,25 @@ int main(int argc, char *argv[])
meshToMesh::interpolationMethod method =
meshToMesh::interpolationMethodNames_[mapMethod];
patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodNames_
[
meshToMesh::interpolationMethodAMI(method)
];
patchMapMethod =
AMIPatchToPatchInterpolation::interpolationMethodNames_
[
meshToMesh::interpolationMethodAMI(method)
];
}
word procMapMethod =
meshToMesh::procMapMethodNames_
[
meshToMesh::procMapMethod::pmAABB
];
if (args.readIfPresent("prcoMapMethod", procMapMethod))
{
Info<< "Processor map method: " << procMapMethod << endl;
}
// Optionally override
if (args.readIfPresent("patchMapMethod", patchMapMethod))
{
@ -325,6 +353,7 @@ int main(int argc, char *argv[])
meshTarget,
mapMethod,
patchMapMethod,
procMapMethod,
subtract,
selectedFields,
noLagrangian
@ -340,6 +369,7 @@ int main(int argc, char *argv[])
cuttingPatches,
mapMethod,
patchMapMethod,
procMapMethod,
subtract,
selectedFields,
noLagrangian

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -65,10 +65,17 @@ void Foam::fv::interRegionOption::setMapper()
(
mesh_,
nbrMesh,
meshToMesh::interpolationMethodNames_.lookup
meshToMesh::interpolationMethodNames_.lookupOrDefault
(
"interpolationMethod",
coeffs_
coeffs_,
meshToMesh::interpolationMethod::imCellVolumeWeight
),
meshToMesh::procMapMethodNames_.lookupOrDefault
(
"procMapMethod",
coeffs_,
meshToMesh::procMapMethod::pmAABB
),
false // not interpolating patches
)

View File

@ -784,9 +784,10 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
(
srcMesh,
tgtMesh,
meshToMesh::imCellVolumeWeight,
meshToMesh::interpolationMethod::imCellVolumeWeight,
HashTable<word>(0), // patchMap,
wordList(0), // cuttingPatches
meshToMesh::procMapMethod::pmAABB,
false // do not normalise
);

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -84,7 +84,14 @@ bool Foam::meshToMeshMethod::intersect
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
// Note: avoid demand-driven construction of cellPoints
// treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
for (label i = 1; i < cellFaces.size(); ++i)
{
bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
}
return overlapEngine.cellCellOverlapMinDecomp
(
@ -106,7 +113,14 @@ Foam::scalar Foam::meshToMeshMethod::interVol
{
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
// Note: avoid demand-driven construction of cellPoints
// treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
for (label i = 1; i < cellFaces.size(); ++i)
{
bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
}
scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
(
@ -124,21 +138,28 @@ Foam::scalar Foam::meshToMeshMethod::interVol
Foam::Tuple2<Foam::scalar, Foam::point>
Foam::meshToMeshMethod::interVolAndCentroid
(
const label srcCellI,
const label tgtCellI
const label srcCelli,
const label tgtCelli
)
{
tetOverlapVolume overlapEngine;
treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCellI]);
// Note: avoid demand-driven construction of cellPoints
// treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
for (label i = 1; i < cellFaces.size(); ++i)
{
bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
}
Tuple2<scalar, point> volAndInertia =
overlapEngine.cellCellOverlapMomentMinDecomp
(
src_,
srcCellI,
srcCelli,
tgt_,
tgtCellI,
tgtCelli,
bbTgtCell
);

View File

@ -52,6 +52,17 @@ Foam::meshToMesh::interpolationMethodNames_
};
const Foam::Enum
<
Foam::meshToMesh::procMapMethod
>
Foam::meshToMesh::procMapMethodNames_
{
{ procMapMethod::pmAABB, "AABB" },
{ procMapMethod::pmLOD, "LOD" },
};
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<>
@ -628,18 +639,18 @@ Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method)
{
switch (method)
{
case imDirect:
case interpolationMethod::imDirect:
{
return AMIPatchToPatchInterpolation::imDirect;
break;
}
case imMapNearest:
case interpolationMethod::imMapNearest:
{
return AMIPatchToPatchInterpolation::imMapNearest;
break;
}
case imCellVolumeWeight:
case imCorrectedCellVolumeWeight:
case interpolationMethod::imCellVolumeWeight:
case interpolationMethod::imCorrectedCellVolumeWeight:
{
return AMIPatchToPatchInterpolation::imFaceAreaWeight;
break;
@ -647,7 +658,7 @@ Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method)
default:
{
FatalErrorInFunction
<< "Unhandled enumeration " << method
<< "Unhandled enumeration " << interpolationMethodNames_[method]
<< abort(FatalError);
}
}
@ -827,11 +838,13 @@ Foam::meshToMesh::meshToMesh
const polyMesh& src,
const polyMesh& tgt,
const interpolationMethod& method,
const procMapMethod& mapMethod,
bool interpAllPatches
)
:
srcRegion_(src),
tgtRegion_(tgt),
procMapMethod_(mapMethod),
srcPatchID_(),
tgtPatchID_(),
patchAMIs_(),
@ -865,11 +878,13 @@ Foam::meshToMesh::meshToMesh
const polyMesh& tgt,
const word& methodName,
const word& AMIMethodName,
const procMapMethod& mapMethod,
bool interpAllPatches
)
:
srcRegion_(src),
tgtRegion_(tgt),
procMapMethod_(mapMethod),
srcPatchID_(),
tgtPatchID_(),
patchAMIs_(),
@ -896,11 +911,13 @@ Foam::meshToMesh::meshToMesh
const interpolationMethod& method,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const procMapMethod& mapMethod,
const bool normalise
)
:
srcRegion_(src),
tgtRegion_(tgt),
procMapMethod_(mapMethod),
srcPatchID_(),
tgtPatchID_(),
patchAMIs_(),
@ -936,11 +953,13 @@ Foam::meshToMesh::meshToMesh
const word& AMIMethodName, // boundary mapping
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const procMapMethod& mapMethod,
const bool normalise
)
:
srcRegion_(src),
tgtRegion_(tgt),
procMapMethod_(mapMethod),
srcPatchID_(),
tgtPatchID_(),
patchAMIs_(),

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -66,7 +66,7 @@ public:
// Public data types
//- Enumeration specifying interpolation method
enum interpolationMethod
enum class interpolationMethod
{
imDirect,
imMapNearest,
@ -76,6 +76,15 @@ public:
static const Enum<interpolationMethod> interpolationMethodNames_;
//- Enumeration specifying processor parallel map construction method
enum class procMapMethod
{
pmAABB,
pmLOD
};
static const Enum<procMapMethod> procMapMethodNames_;
private:
// Private data
@ -86,6 +95,9 @@ private:
//- Reference to the target mesh
const polyMesh& tgtRegion_;
//- Processor parallel map construction method
procMapMethod procMapMethod_;
//- List of target patch IDs per source patch (local index)
List<label> srcPatchID_;
@ -301,6 +313,7 @@ public:
const polyMesh& src,
const polyMesh& tgt,
const interpolationMethod& method,
const procMapMethod& mapMethod = procMapMethod::pmAABB,
const bool interpAllPatches = true
);
@ -311,6 +324,7 @@ public:
const polyMesh& tgt,
const word& methodName, // internal mapping
const word& AMIMethodName, // boundary mapping
const procMapMethod& mapMethod = procMapMethod::pmAABB,
const bool interpAllPatches = true
);
@ -322,6 +336,7 @@ public:
const interpolationMethod& method,
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const procMapMethod& mapMethod = procMapMethod::pmAABB,
const bool normalise = true
);
@ -335,6 +350,7 @@ public:
const word& AMIMethodName, // boundary mapping
const HashTable<word>& patchMap,
const wordList& cuttingPatches,
const procMapMethod& mapMethod = procMapMethod::pmAABB,
const bool normalise = true
);

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -31,6 +31,7 @@ License
#include "processorPolyPatch.H"
#include "SubField.H"
#include "AABBTree.H"
#include "cellBox.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -98,12 +99,12 @@ Foam::label Foam::meshToMesh::calcOverlappingProcs
{
const treeBoundBoxList& bbp = procBb[proci];
forAll(bbp, bbi)
for (const treeBoundBox& b : bbp)
{
if (bbp[bbi].overlaps(bb))
if (b.overlaps(bb))
{
overlaps[proci] = true;
nOverlaps++;
++nOverlaps;
break;
}
}
@ -119,139 +120,179 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
const polyMesh& tgt
) const
{
// get decomposition of cells on src mesh
List<treeBoundBoxList> procBb(Pstream::nProcs());
if (src.nCells() > 0)
switch (procMapMethod_)
{
procBb[Pstream::myProcNo()] = AABBTree<labelList>
(
src.cellPoints(),
src.points(),
false
).boundBoxes();
}
else
{
procBb[Pstream::myProcNo()] = treeBoundBoxList();
}
Pstream::gatherList(procBb);
Pstream::scatterList(procBb);
if (debug)
{
InfoInFunction
<< "Determining extent of src mesh per processor:" << nl
<< "\tproc\tbb" << endl;
forAll(procBb, proci)
case procMapMethod::pmLOD:
{
Info<< '\t' << proci << '\t' << procBb[proci] << endl;
Info<< "meshToMesh: Using processorLOD method" << endl;
// Create processor map of overlapping faces. This map gets
// (possibly remote) cells from the tgt mesh such that they
// (together) cover all of the src mesh
label nGlobalSrcCells =
returnReduce(src.cells().size(), sumOp<label>());
label cellsPerBox = max(1, 0.001*nGlobalSrcCells);
typename processorLODs::cellBox boxLOD
(
src.cells(),
src.faces(),
src.points(),
tgt.cells(),
tgt.faces(),
tgt.points(),
cellsPerBox,
src.nCells()
);
return boxLOD.map();
break;
}
}
// determine which cells of tgt mesh overlaps src mesh per proc
const cellList& cells = tgt.cells();
const faceList& faces = tgt.faces();
const pointField& points = tgt.points();
labelListList sendMap;
{
// per processor indices into all segments to send
List<DynamicList<label>> dynSendMap(Pstream::nProcs());
label iniSize = floor(tgt.nCells()/Pstream::nProcs());
forAll(dynSendMap, proci)
default:
{
dynSendMap[proci].setCapacity(iniSize);
}
Info<< "meshToMesh: Using AABBTree method" << endl;
// work array - whether src processor bb overlaps the tgt cell bounds
boolList procBbOverlaps(Pstream::nProcs());
forAll(cells, celli)
{
const cell& c = cells[celli];
// get decomposition of cells on src mesh
List<treeBoundBoxList> procBb(Pstream::nProcs());
// determine bounding box of tgt cell
boundBox cellBb(boundBox::invertedBox);
forAll(c, facei)
if (src.nCells() > 0)
{
const face& f = faces[c[facei]];
forAll(f, fp)
procBb[Pstream::myProcNo()] = AABBTree<labelList>
(
src.cellPoints(),
src.points(),
false
).boundBoxes();
}
else
{
procBb[Pstream::myProcNo()] = treeBoundBoxList();
}
Pstream::gatherList(procBb);
Pstream::scatterList(procBb);
if (debug)
{
InfoInFunction
<< "Determining extent of src mesh per processor:" << nl
<< "\tproc\tbb" << endl;
forAll(procBb, proci)
{
cellBb.add(points, f);
Info<< '\t' << proci << '\t' << procBb[proci] << endl;
}
}
// find the overlapping tgt cells on each src processor
(void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
forAll(procBbOverlaps, proci)
// determine which cells of tgt mesh overlaps src mesh per proc
const cellList& cells = tgt.cells();
const faceList& faces = tgt.faces();
const pointField& points = tgt.points();
labelListList sendMap;
{
if (procBbOverlaps[proci])
// per processor indices into all segments to send
List<DynamicList<label>> dynSendMap(Pstream::nProcs());
label iniSize = floor(tgt.nCells()/Pstream::nProcs());
forAll(dynSendMap, proci)
{
dynSendMap[proci].append(celli);
dynSendMap[proci].setCapacity(iniSize);
}
// work array - whether src processor bb overlaps the tgt cell
// bounds
boolList procBbOverlaps(Pstream::nProcs());
forAll(cells, celli)
{
const cell& c = cells[celli];
// determine bounding box of tgt cell
boundBox cellBb(boundBox::invertedBox);
forAll(c, facei)
{
const face& f = faces[c[facei]];
forAll(f, fp)
{
cellBb.add(points, f);
}
}
// find the overlapping tgt cells on each src processor
(void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
forAll(procBbOverlaps, proci)
{
if (procBbOverlaps[proci])
{
dynSendMap[proci].append(celli);
}
}
}
// convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
sendMap[proci].transfer(dynSendMap[proci]);
}
}
}
// convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
sendMap[proci].transfer(dynSendMap[proci]);
// debug printing
if (debug)
{
Pout<< "Of my " << cells.size()
<< " target cells I need to send to:" << nl
<< "\tproc\tcells" << endl;
forAll(sendMap, proci)
{
Pout<< '\t' << proci << '\t' << sendMap[proci].size()
<< endl;
}
}
// send over how many tgt cells I need to receive from each
// processor
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
}
Pstream::gatherList(sendSizes);
Pstream::scatterList(sendSizes);
// determine order of receiving
labelListList constructMap(Pstream::nProcs());
label segmentI = 0;
forAll(constructMap, proci)
{
// what I need to receive is what other processor is sending
// to me
label nRecv = sendSizes[proci][Pstream::myProcNo()];
constructMap[proci].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[proci][i] = segmentI++;
}
}
return autoPtr<mapDistribute>::New
(
segmentI, // size after construction
std::move(sendMap),
std::move(constructMap)
);
break;
}
}
// debug printing
if (debug)
{
Pout<< "Of my " << cells.size() << " target cells I need to send to:"
<< nl << "\tproc\tcells" << endl;
forAll(sendMap, proci)
{
Pout<< '\t' << proci << '\t' << sendMap[proci].size() << endl;
}
}
// send over how many tgt cells I need to receive from each processor
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
}
Pstream::gatherList(sendSizes);
Pstream::scatterList(sendSizes);
// determine order of receiving
labelListList constructMap(Pstream::nProcs());
label segmentI = 0;
forAll(constructMap, proci)
{
// what I need to receive is what other processor is sending to me
label nRecv = sendSizes[proci][Pstream::myProcNo()];
constructMap[proci].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[proci][i] = segmentI++;
}
}
return autoPtr<mapDistribute>::New
(
segmentI, // size after construction
std::move(sendMap),
std::move(constructMap)
);
}