ENH: Initial parallelisation - point exchange over bound box interfaces.

This commit is contained in:
graham
2011-03-10 13:42:41 +00:00
parent 8c3bb1eee6
commit eede286417
11 changed files with 405 additions and 60 deletions

View File

@ -514,7 +514,7 @@ void Foam::conformalVoronoiMesh::createFeaturePoints()
insertMixedFeaturePoints();
Info<< " Inserted " << number_of_vertices() << " vertices" << endl;
Pout<< " Inserted " << number_of_vertices() << " vertices" << endl;
featureVertices_.setSize(number_of_vertices());
@ -563,10 +563,15 @@ void Foam::conformalVoronoiMesh::insertConvexFeaturePoints()
ptI++
)
{
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
const Foam::point& featPt = feMesh.points()[ptI];
if (!geometryToConformTo_.positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
@ -610,10 +615,15 @@ void Foam::conformalVoronoiMesh::insertConcaveFeaturePoints()
ptI++
)
{
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
const Foam::point& featPt = feMesh.points()[ptI];
if (!geometryToConformTo_.positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
@ -703,6 +713,11 @@ void Foam::conformalVoronoiMesh::insertMixedFeaturePoints()
const Foam::point& pt(feMesh.points()[ptI]);
if (!geometryToConformTo_.positionOnThisProc(pt))
{
continue;
}
scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
forAll(pEds, e)
@ -1162,7 +1177,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
:
Delaunay(),
runTime_(runTime),
rndGen_(7864293),
rndGen_(64293*Pstream::myProcNo()),
allGeometry_
(
IOobject

View File

@ -416,6 +416,18 @@ private:
const Delaunay::Finite_vertices_iterator& vit
) const;
//- Check to see if dual cell specified by given vertex iterator
// intersects the parallel interface, hence needs referred. Returns
// destination processor(s).
List<label> parallelInterfaceIntersection
(
const Delaunay::Finite_vertices_iterator& vit
) const;
//- Which processor should a vertex generating the given dual edge be
// referred to
label targetProc(direction ptOnFace) const;
//- Find the "worst" protrusion of a dual cell through the surface,
// subject to the maxSurfaceProtrusion tolerance
void dualCellLargestSurfaceProtrusion

View File

@ -152,6 +152,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
hitSurface
);
// Pout<< "vert " << vert << endl;
if (surfHit.hit())
{
vit->setNearBoundary();
@ -176,10 +178,12 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
);
}
}
// Pout<< " surfHit " << surfHit << endl;
}
}
Info<< nl << "Initial conformation" << nl
Pout<< nl << "Initial conformation" << nl
<< " Number of vertices " << number_of_vertices() << nl
<< " Number of surface hits " << surfaceHits.size() << nl
<< " Number of edge hits " << featureEdgeHits.size()
@ -204,6 +208,171 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
initialTotalHits = surfaceHits.size() + featureEdgeHits.size();
}
label totalInterfacePoints = 0;
label nIter = 0;
if (Pstream::parRun())
{
do
{
// Propagating vertices to other processors to form halo information
DynamicList<Foam::point> parallelInterfacePoints;
DynamicList<label> targetProcessor;
for
(
Delaunay::Finite_vertices_iterator vit =
finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->internalOrBoundaryPoint())
{
List<label> toProcs = parallelInterfaceIntersection(vit);
forAll(toProcs, tPI)
{
label toProc = toProcs[tPI];
if (toProc > -1)
{
parallelInterfacePoints.append
(
topoint(vit->point())
);
targetProcessor.append(toProc);
std::list<Vertex_handle> incidentVertices;
incident_vertices
(
vit,
std::back_inserter(incidentVertices)
);
for
(
std::list<Vertex_handle>::iterator ivit =
incidentVertices.begin();
ivit != incidentVertices.end();
++ivit
)
{
parallelInterfacePoints.append
(
topoint((*ivit)->point())
);
targetProcessor.append(toProc);
}
}
}
}
}
writePoints("parallelInterfacePoints.obj", parallelInterfacePoints);
// Determine send map
// ~~~~~~~~~~~~~~~~~~
// 1. Count
labelList nSend(Pstream::nProcs(), 0);
forAll(targetProcessor, i)
{
label procI = targetProcessor[i];
nSend[procI]++;
}
// Send over how many I need to receive
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()] = nSend;
combineReduce(sendSizes, UPstream::listEq());
// 2. Size sendMap
labelListList sendMap(Pstream::nProcs());
forAll(nSend, procI)
{
sendMap[procI].setSize(nSend[procI]);
nSend[procI] = 0;
}
// 3. Fill sendMap
forAll(targetProcessor, i)
{
label procI = targetProcessor[i];
sendMap[procI][nSend[procI]++] = i;
}
// Determine receive map
// ~~~~~~~~~~~~~~~~~~~~~
labelListList constructMap(Pstream::nProcs());
// Local transfers first
constructMap[Pstream::myProcNo()] = identity
(
sendMap[Pstream::myProcNo()].size()
);
label constructSize = constructMap[Pstream::myProcNo()].size();
forAll(constructMap, procI)
{
if (procI != Pstream::myProcNo())
{
label nRecv = sendSizes[procI][Pstream::myProcNo()];
constructMap[procI].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[procI][i] = constructSize++;
}
}
}
mapDistribute pointMap
(
constructSize,
sendMap.xfer(),
constructMap.xfer()
);
totalInterfacePoints = parallelInterfacePoints.size();
reduce(totalInterfacePoints, sumOp<label>());
pointMap.distribute(parallelInterfacePoints);
forAll(parallelInterfacePoints, ptI)
{
insertPoint
(
parallelInterfacePoints[ptI],
Vb::ptFarPoint
);
}
Info<< "totalInterfacePoints " << totalInterfacePoints << endl;
nIter++;
} while (totalInterfacePoints > 0 && nIter < 3);
}
label iterationNo = 0;
label maxIterations =
@ -214,7 +383,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
label hitLimit = label(iterationToIntialHitRatioLimit*initialTotalHits);
Info<< nl << "Stopping iterations when: " << nl
Pout<< nl << "Stopping iterations when: " << nl
<<" total number of hits drops below "
<< iterationToIntialHitRatioLimit << " of initial hits ("
<< hitLimit << ")" << nl
@ -469,7 +638,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
// }
// }
storeSurfaceConformation();
Pout<< "NOT STORING SURFACE CONFORMATION" << endl;
//storeSurfaceConformation();
}
@ -498,21 +668,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
}
Foam::point dE0 = topoint(dual(fit->first));
// If edge end is outside bounding box then edge cuts boundary
if (!geometryToConformTo_.bounds().contains(dE0))
{
return true;
}
Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
// If other edge end is outside bounding box then edge cuts boundary
if (!geometryToConformTo_.bounds().contains(dE1))
{
return true;
}
// Check for the edge passing through a surface
if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
{
@ -524,6 +681,115 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
}
Foam::List<Foam::label>
Foam::conformalVoronoiMesh::parallelInterfaceIntersection
(
const Delaunay::Finite_vertices_iterator& vit
) const
{
std::list<Facet> facets;
incident_facets(vit, std::back_inserter(facets));
DynamicList<label> procs;
for
(
std::list<Facet>::iterator fit=facets.begin();
fit != facets.end();
++fit
)
{
if
(
is_infinite(fit->first)
|| is_infinite(fit->first->neighbor(fit->second))
)
{
return procs;
}
Foam::point dE0 = topoint(dual(fit->first));
Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
Foam::point boxPt(vector::one*GREAT);
direction ptOnFace = -1;
geometryToConformTo_.bounds().intersects
(
dE0,
dE1 - dE0,
dE0,
dE1,
boxPt,
ptOnFace
);
if
(
ptOnFace > 0
&& !geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1)
)
{
// A surface penetration is not found, and a bounds penetration is
// found, so the dual edge has crossed the parallel interface.
label target = targetProc(ptOnFace);
if
(
target > -1
&& findIndex(procs, target) == -1
)
{
procs.append(target);
}
// Pout<< "Parallel box penetration, face " << ptOnFace
// << " proc " << target << endl;
// meshTools::writeOBJ(Pout, dE0);
// meshTools::writeOBJ(Pout, dE1);
// meshTools::writeOBJ(Pout, boxPt);
// Pout << "l 1 2" << endl;
}
}
return procs;
}
Foam::label Foam::conformalVoronoiMesh::targetProc
(
direction ptOnFace
) const
{
// Hard coded to 8 proc, bound-box octant split:
label faceIndex = 0;
// From: http://graphics.stanford.edu/~seander/bithacks.html
while (ptOnFace >>= 1)
{
faceIndex++;
}
label procArray[8][6] =
{
{-1, 1, -1, 2, -1, 4},
{ 0, -1, -1, 3, -1, 5},
{-1, 3, 0, -1, -1, 6},
{ 2, -1, 1, -1, -1, 7},
{-1, 5, -1, 6, 0, -1},
{ 4, -1, -1, 7, 1, -1},
{-1, 7, 4, -1, 2, -1},
{ 6, -1, 5, -1, 3, -1}
};
return initListList<List<label>, label, 8, 6>(procArray)
[Pstream::myProcNo()][faceIndex];
}
void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
(
const Delaunay::Finite_vertices_iterator& vit,

View File

@ -94,6 +94,11 @@ bool Foam::conformalVoronoiMesh::insertSpecialisedFeaturePoint
if (nExternal == 2 && nInternal == 1)
{
// if (!geometryToConformTo_.positionOnThisProc(pt))
// {
// continue;
// }
// Info<< "nExternal == 2 && nInternal == 1" << endl;
// const Foam::point& featPt = feMesh.points()[ptI];

View File

@ -240,14 +240,14 @@ inline void Foam::conformalVoronoiMesh::insertPoints
const std::vector<Point>& points
)
{
Info<< " " << points.size() << " points to insert..." << endl;
Pout<< " " << points.size() << " points to insert..." << endl;
label nVert = number_of_vertices();
// using the range insert (faster than inserting points one by one)
insert(points.begin(), points.end());
Info<< " " << number_of_vertices() - startOfInternalPoints_
Pout<< " " << number_of_vertices() - startOfInternalPoints_
<< " points inserted" << endl;
for

View File

@ -150,7 +150,7 @@ void Foam::conformalVoronoiMesh::writeInternalDelaunayVertices
internalDelaunayVertices
);
Info<< nl
Pout<< nl
<< "Writing " << internalDVs.name()
<< " to " << internalDVs.instance()
<< endl;
@ -228,7 +228,7 @@ void Foam::conformalVoronoiMesh::writeMesh
filterFaces
);
Info<< nl << "Writing polyMesh to " << instance << endl;
Pout<< nl << "Writing polyMesh to " << instance << endl;
writeMesh
(
@ -572,7 +572,7 @@ void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
if (!protrudingCells.empty())
{
Info<< nl << "Found " << protrudingCells.size()
Pout<< nl << "Found " << protrudingCells.size()
<< " cells protruding from the surface, writing cellSet "
<< protrudingCells.name()
<< endl;

View File

@ -59,6 +59,7 @@ class indexedVertex
// ptNearBoundaryPoint : internal near boundary point.
// ptInternalPoint : internal point.
// ptFarPoint : far-point.
// > ptFarPoint, < 0 : referred point from processor -(type_ + 1)
// >= 0 : part of point-pair. Index of other point.
// Lowest numbered is inside one (master).
int type_;
@ -72,10 +73,9 @@ public:
enum pointTypes
{
ptNearBoundaryPoint = -4,
ptInternalPoint = -3,
ptMirrorPoint = -2,
ptFarPoint = -1
ptNearBoundaryPoint = INT_MIN,
ptInternalPoint = INT_MIN + 1,
ptFarPoint = INT_MIN + 2
};
typedef typename Vb::Vertex_handle Vertex_handle;
@ -198,6 +198,27 @@ public:
}
// is this a referred vertex
inline bool referred() const
{
return type_ < 0 && type_ > ptFarPoint;
}
// For referred vertices, what is the processor index
inline int procIndex() const
{
if (referred())
{
return -(type_ + 1);
}
else
{
return -1;
}
}
//- Set the point to be internal
inline void setInternal()
{
@ -218,12 +239,6 @@ public:
type_ = ptNearBoundaryPoint;
}
//- Is point a mirror point
inline bool mirrorPoint() const
{
return type_ == ptMirrorPoint;
}
//- Either master or slave of pointPair.
inline bool pairPoint() const
@ -270,7 +285,7 @@ public:
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const
{
return pairPoint() || mirrorPoint() || nearBoundary();
return pairPoint() || nearBoundary();
}

View File

@ -361,6 +361,15 @@ bool Foam::conformationSurfaces::overlaps(const treeBoundBox& bb) const
}
bool Foam::conformationSurfaces::positionOnThisProc(const point& pt) const
{
// This is likely to give problems when a point is on the boundary between
// two processors.
return bounds_.contains(pt);
}
Foam::Field<bool> Foam::conformationSurfaces::inside
(
const pointField& samplePts

View File

@ -170,6 +170,9 @@ public:
// the surfaces
bool overlaps(const treeBoundBox& bb) const;
//- Check if the point is in the domain handled by this processor
bool positionOnThisProc(const point& pt) const;
//- Check if points are inside surfaces to conform to
Field<bool> inside(const pointField& samplePts) const;

View File

@ -119,7 +119,7 @@ void Foam::hierarchicalDensityWeightedStochastic::recurseAndFill
if (debug)
{
Info<< newName + "_overlap " << subBB << endl;
Pout<< newName + "_overlap " << subBB << endl;
}
if (!fillBox(initialPoints, subBB, true))
@ -144,7 +144,7 @@ void Foam::hierarchicalDensityWeightedStochastic::recurseAndFill
if (debug)
{
Info<< newName + "_inside " << subBB << endl;
Pout<< newName + "_inside " << subBB << endl;
}
if (!fillBox(initialPoints, subBB, false))
@ -221,7 +221,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (surfHit.hit())
{
// Info<< "box wellInside, no need to sample surface." << endl;
// Pout<< "box wellInside, no need to sample surface." << endl;
wellInside = true;
}
@ -248,7 +248,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
minimumSurfaceDistanceCoeffSqr_*sqr(cornerSizes)
);
// Info<< corners << nl << cornerSizes << nl << insideCorners << endl;
// Pout<< corners << nl << cornerSizes << nl << insideCorners << endl;
forAll(insideCorners, i)
{
@ -267,7 +267,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (maxCellSize/minCellSize > maxSizeRatio_)
{
// Info<< "Abort fill at corner sample stage,"
// Pout<< "Abort fill at corner sample stage,"
// << " minCellSize " << minCellSize
// << " maxCellSize " << maxCellSize
// << " maxSizeRatio " << maxCellSize/minCellSize
@ -281,7 +281,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
// If one or more corners is not "wellInside", then treat this
// as an overlapping box.
// Info<< "Inside box found to have some non-wellInside "
// Pout<< "Inside box found to have some non-wellInside "
// << "corners, using overlapping fill."
// << endl;
@ -378,7 +378,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (maxCellSize/minCellSize > maxSizeRatio_)
{
// Info<< "Abort fill at surface sample stage, "
// Pout<< "Abort fill at surface sample stage, "
// << " minCellSize " << minCellSize
// << " maxCellSize " << maxCellSize
// << " maxSizeRatio " << maxCellSize/minCellSize
@ -393,7 +393,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
// then treat this as an overlapping box.
overlapping = true;
// Info<< "Inside box found to have some non-"
// Pout<< "Inside box found to have some non-"
// << "wellInside surface points, using "
// << "overlapping fill."
// << endl;
@ -480,7 +480,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (maxCellSize/minCellSize > maxSizeRatio_)
{
// Info<< "Abort fill at sample stage,"
// Pout<< "Abort fill at sample stage,"
// << " minCellSize " << minCellSize
// << " maxCellSize " << maxCellSize
// << " maxSizeRatio " << maxCellSize/minCellSize
@ -493,17 +493,17 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (nInside == 0)
{
// Info<< "No sample points found inside box" << endl;
// Pout<< "No sample points found inside box" << endl;
return true;
}
// Info<< scalar(nInside)/scalar(samplePoints.size())
// Pout<< scalar(nInside)/scalar(samplePoints.size())
// << " full overlapping box" << endl;
totalVolume *= scalar(nInside)/scalar(samplePoints.size());
// Info<< "Total volume to fill = " << totalVolume << endl;
// Pout<< "Total volume to fill = " << totalVolume << endl;
// Using the sampledPoints as the first test locations as they are
// randomly shuffled, but unfiormly sampling space and have wellInside
@ -547,7 +547,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
scalar r = rnd.scalar01();
// Info<< "totalVolume " << totalVolume << nl
// Pout<< "totalVolume " << totalVolume << nl
// << "volumeAdded " << volumeAdded << nl
// << "localVolume " << localVolume << nl
// << "addProbability " << addProbability << nl
@ -559,7 +559,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
// Place this volume before finishing filling this
// box
// Info<< "Final volume probability break accept"
// Pout<< "Final volume probability break accept"
// << endl;
initialPoints.push_back
@ -583,7 +583,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (volumeAdded < totalVolume)
{
// Info<< "Adding random points, remaining volume "
// Pout<< "Adding random points, remaining volume "
// << totalVolume - volumeAdded
// << endl;
@ -633,7 +633,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
if (maxCellSize/minCellSize > maxSizeRatio_)
{
// Info<< "Abort fill at random fill stage,"
// Pout<< "Abort fill at random fill stage,"
// << " minCellSize " << minCellSize
// << " maxCellSize " << maxCellSize
// << " maxSizeRatio " << maxCellSize/minCellSize
@ -665,7 +665,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
scalar r = rnd.scalar01();
// Info<< "totalVolume " << totalVolume << nl
// Pout<< "totalVolume " << totalVolume << nl
// << "volumeAdded " << volumeAdded << nl
// << "localVolume " << localVolume << nl
// << "addProbability " << addProbability << nl
@ -677,7 +677,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
// Place this volume before finishing filling this
// box
// Info<< "Final volume probability break accept"
// Pout<< "Final volume probability break accept"
// << endl;
initialPoints.push_back
@ -701,7 +701,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
globalTrialPoints_ += trialPoints;
// Info<< trialPoints
// Pout<< trialPoints
// << " locations queried, " << initialPoints.size() - initialSize
// << " points placed, ("
// << scalar(initialPoints.size() - initialSize)
@ -777,7 +777,7 @@ hierarchicalDensityWeightedStochastic::initialPoints() const
if (debug)
{
Info<< " Filling box " << hierBB << endl;
Pout<< " Filling box " << hierBB << endl;
}
recurseAndFill
@ -788,10 +788,18 @@ hierarchicalDensityWeightedStochastic::initialPoints() const
"recursionBox"
);
Info<< " " << initialPoints.size() << " points placed" << nl
label nInitialPoints = initialPoints.size();
if (Pstream::parRun())
{
reduce(nInitialPoints, sumOp<label>());
reduce(globalTrialPoints_, sumOp<label>());
}
Info<< " " << nInitialPoints << " points placed" << nl
<< " " << globalTrialPoints_ << " locations queried" << nl
<< " "
<< scalar(initialPoints.size())/scalar(max(globalTrialPoints_, 1))
<< scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
<< " success rate"
<< endl;

View File

@ -74,6 +74,18 @@ std::vector<Vb::Point> pointFile::initialPoints() const
<< exit(FatalError) << endl;
}
boolList procPt(points.size(), false);
forAll(points, ptI)
{
procPt[ptI] = cvMesh_.geometryToConformTo().positionOnThisProc
(
points[ptI]
);
}
inplaceSubset(procPt, points);
std::vector<Vb::Point> initialPoints;
Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside