ENH: Parallel face collapse WIP.

This commit is contained in:
graham
2011-06-30 18:42:23 +01:00
parent f7f069083d
commit 86514c689c
5 changed files with 363 additions and 120 deletions

View File

@ -83,6 +83,11 @@ class faceAreaWeightModel;
class backgroundMeshDecomposition;
typedef FixedList<label, 4>globalTetIndices;
typedef HashTable<label, globalTetIndices, Hash<globalTetIndices> >
HashTableGlobalTetIndices;
/*---------------------------------------------------------------------------*\
Class conformalVoronoiMesh Declaration
\*---------------------------------------------------------------------------*/
@ -702,6 +707,13 @@ private:
const Delaunay::Finite_edges_iterator& eit
) const;
//- Which processors are attached to the dual edge represented by this
// Delaunay facet
inline List<label> processorsAttached
(
const Delaunay::Finite_facets_iterator& fit
) const;
//- Merge vertices that are very close together
void mergeCloseDualVertices
(
@ -783,8 +795,7 @@ private:
void indexDualVertices
(
pointField& pts,
PackedBoolList& boundaryPts,
const globalIndex& globalParallelPointIndices
PackedBoolList& boundaryPts
);
//- Re-index all of the the Delaunay cells

View File

@ -25,6 +25,7 @@ License
#include "conformalVoronoiMesh.H"
#include "motionSmoother.H"
#include "backgroundMeshDecomposition.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -76,7 +77,8 @@ void Foam::conformalVoronoiMesh::calcDualMesh
// || !cit->vertex(3)->real()
// )
// {
// cit->filterCount() = 100;
// cit->filterCount() =
// cvMeshControls().filterCountSkipThreshold() + 1;
// }
// else
// {
@ -126,7 +128,8 @@ void Foam::conformalVoronoiMesh::calcDualMesh
// ++cit
// )
// {
// (*cit)->filterCount() = 100;
// (*cit)->filterCount() =
// cvMeshControls().filterCountSkipThreshold() + 1;
// }
// }
// }
@ -134,13 +137,10 @@ void Foam::conformalVoronoiMesh::calcDualMesh
PackedBoolList boundaryPts(number_of_cells(), false);
globalIndex globalParallelPointIndices(number_of_cells());
indexDualVertices
(
points,
boundaryPts,
globalParallelPointIndices
boundaryPts
);
{
@ -202,8 +202,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
indexDualVertices
(
points,
boundaryPts,
globalParallelPointIndices
boundaryPts
);
{
@ -509,7 +508,12 @@ void Foam::conformalVoronoiMesh::mergeCloseDualVertices
{
Map<label> dualPtIndexMap;
nPtsMerged = mergeCloseDualVertices(pts, boundaryPts, dualPtIndexMap);
nPtsMerged = mergeCloseDualVertices
(
pts,
boundaryPts,
dualPtIndexMap
);
// if (nPtsMerged > 0)
// {
@ -522,9 +526,12 @@ void Foam::conformalVoronoiMesh::mergeCloseDualVertices
} while (nPtsMerged > 0);
Info<< " Merged "
<< returnReduce(nPtsMergedSum, sumOp<label>())
<< " points " << endl;
reduce(nPtsMergedSum, sumOp<label>());
if (nPtsMergedSum > 0)
{
Info<< " Merged " << nPtsMergedSum << " points " << endl;
}
}
@ -539,6 +546,35 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
scalar closenessTolerance = cvMeshControls().mergeClosenessCoeff();
globalIndex globalDelaunayVertexIndices(number_of_vertices());
HashTableGlobalTetIndices dualVertexGlobalTetIndices;
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
cit != finite_cells_end();
++cit
)
{
if
(
cit->internalOrBoundaryDualVertex()
&& cit->parallelDualVertex()
)
{
dualVertexGlobalTetIndices.insert
(
cit->vertexGlobalIndices(globalDelaunayVertexIndices),
cit->cellIndex()
);
}
}
DynamicList<label> toProc;
DynamicList<globalTetIndices> globalToMapFirst;
DynamicList<globalTetIndices> globalToMapSecond;
for
(
Delaunay::Finite_facets_iterator fit = finite_facets_begin();
@ -569,26 +605,134 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
< sqr(averageAnyCellSize(fit)*closenessTolerance)
)
{
if (boundaryPts[c2I] == true)
if (c1->parallelDualVertex() && c2->parallelDualVertex())
{
// Both dual vertices are on a processor patch face
// The collapse is actioned by the master processor
label masterProc = min
(
c1->dualVertexMasterProc(),
c2->dualVertexMasterProc()
);
if (masterProc == Pstream::myProcNo())
{
// perform the collapse and insert values to notify
// everyone else
labelList attachedProcs(processorsAttached(fit));
forAll(attachedProcs, aPI)
{
toProc.append(attachedProcs[aPI]);
if (boundaryPts[c2I] == true)
{
// If c2I is a boundary point, then it is kept.
// If both are boundary points then c2I is
// chosen arbitrarily to be kept.
globalToMapFirst.append
(
c2->vertexGlobalIndices
(
globalDelaunayVertexIndices
)
);
globalToMapSecond.append
(
c1->vertexGlobalIndices
(
globalDelaunayVertexIndices
)
);
}
else
{
globalToMapFirst.append
(
c1->vertexGlobalIndices
(
globalDelaunayVertexIndices
)
);
globalToMapSecond.append
(
c2->vertexGlobalIndices
(
globalDelaunayVertexIndices
)
);
}
}
nPtsMerged++;
}
}
else if (c1->parallelDualVertex() || c2->parallelDualVertex())
{
// If c2I is a boundary point, then it is kept.
// If both are boundary points then c2I is chosen
// arbitrarily to be kept.
dualPtIndexMap.insert(c1I, c2I);
dualPtIndexMap.insert(c2I, c2I);
}
else
{
dualPtIndexMap.insert(c1I, c1I);
dualPtIndexMap.insert(c2I, c1I);
}
nPtsMerged++;
// nPtsMerged++;
}
}
}
mapDistribute map = backgroundMeshDecomposition::buildMap(toProc);
map.distribute(globalToMapFirst);
map.distribute(globalToMapSecond);
forAll(globalToMapFirst, gTMI)
{
label c1I = dualVertexGlobalTetIndices.find(globalToMapFirst[gTMI])();
label c2I = dualVertexGlobalTetIndices.find(globalToMapSecond[gTMI])();
Pout<< globalToMapFirst[gTMI] << " " << c1I << endl;
Pout<< globalToMapSecond[gTMI] << " " << c2I << endl;
dualPtIndexMap.insert(c1I, c1I);
dualPtIndexMap.insert(c2I, c1I);
Pout<< dualPtIndexMap << endl;
}
// if (!c1->farCell() && !c2->farCell() && (c1I != c2I))
// {
// if
// (
// magSqr(pts[c1I] - pts[c2I])
// < sqr(averageAnyCellSize(fit)*closenessTolerance)
// )
// {
// if (boundaryPts[c2I] == true)
// {
// // If c2I is a boundary point, then it is kept.
// // If both are boundary points then c2I is chosen
// // arbitrarily to be kept.
// dualPtIndexMap.insert(c1I, c2I);
// dualPtIndexMap.insert(c2I, c2I);
// }
// else
// {
// dualPtIndexMap.insert(c1I, c1I);
// dualPtIndexMap.insert(c2I, c1I);
// }
// nPtsMerged++;
// }
// }
return nPtsMerged;
}
@ -1551,65 +1695,65 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
List<polyPatch*> patches(patchStarts.size());
// label nValidPatches = 0;
// forAll(patches, p)
// {
// if (patchTypes[p] == processorPolyPatch::typeName)
// {
// // Do not create empty processor patches
// if (patchSizes[p] > 0)
// {
// patches[nValidPatches] = new processorPolyPatch
// (
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// nValidPatches,
// pMesh.boundaryMesh(),
// Pstream::myProcNo(),
// procNeighbours[p]
// );
// nValidPatches++;
// }
// }
// else
// {
// patches[nValidPatches] = polyPatch::New
// (
// patchTypes[p],
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// nValidPatches,
// pMesh.boundaryMesh()
// ).ptr();
// nValidPatches++;
// }
// }
// patches.setSize(nValidPatches);
// pMesh.addPatches(patches);
Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
label nValidPatches = 0;
forAll(patches, p)
{
patches[p] = new polyPatch
(
patchNames[p],
patchSizes[p],
patchStarts[p],
p,
pMesh.boundaryMesh()
);
if (patchTypes[p] == processorPolyPatch::typeName)
{
// Do not create empty processor patches
if (patchSizes[p] > 0)
{
patches[nValidPatches] = new processorPolyPatch
(
patchNames[p],
patchSizes[p],
patchStarts[p],
nValidPatches,
pMesh.boundaryMesh(),
Pstream::myProcNo(),
procNeighbours[p]
);
nValidPatches++;
}
}
else
{
patches[nValidPatches] = polyPatch::New
(
patchTypes[p],
patchNames[p],
patchSizes[p],
patchStarts[p],
nValidPatches,
pMesh.boundaryMesh()
).ptr();
nValidPatches++;
}
}
pMesh.addPatches(patches, false);
patches.setSize(nValidPatches);
pMesh.addPatches(patches);
// Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
// forAll(patches, p)
// {
// patches[p] = new polyPatch
// (
// patchNames[p],
// patchSizes[p],
// patchStarts[p],
// p,
// pMesh.boundaryMesh()
// );
// }
// pMesh.addPatches(patches, false);
// pMesh.overrideCellCentres(cellCentres);
@ -1807,8 +1951,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
void Foam::conformalVoronoiMesh::indexDualVertices
(
pointField& pts,
PackedBoolList& boundaryPts,
const globalIndex& globalParallelPointIndices
PackedBoolList& boundaryPts
)
{
// Indexing Delaunay cells, which are the dual vertices
@ -1821,16 +1964,6 @@ void Foam::conformalVoronoiMesh::indexDualVertices
boundaryPts = false;
OFstream tetStr
(
runTime_.path()
/"processor_"
+ name(Pstream::myProcNo())
+ "_parallelTets.obj"
);
label pTetI = 0;
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@ -1840,35 +1973,23 @@ void Foam::conformalVoronoiMesh::indexDualVertices
{
if (cit->internalOrBoundaryDualVertex())
{
if (cit->parallelDualVertex())
{
cit->cellIndex() = globalParallelPointIndices.toGlobal
(
dualVertI
);
drawDelaunayCell(tetStr, cit, pTetI++);
}
else
{
cit->cellIndex() = dualVertI;
if
(
!cit->vertex(0)->internalOrBoundaryPoint()
|| !cit->vertex(1)->internalOrBoundaryPoint()
|| !cit->vertex(2)->internalOrBoundaryPoint()
|| !cit->vertex(3)->internalOrBoundaryPoint()
)
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
}
}
cit->cellIndex() = dualVertI;
pts[dualVertI] = topoint(dual(cit));
if
(
!cit->vertex(0)->internalOrBoundaryPoint()
|| !cit->vertex(1)->internalOrBoundaryPoint()
|| !cit->vertex(2)->internalOrBoundaryPoint()
|| !cit->vertex(3)->internalOrBoundaryPoint()
)
{
// This is a boundary dual vertex
boundaryPts[dualVertI] = true;
}
dualVertI++;
}
else

View File

@ -347,6 +347,38 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
}
inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
(
const Delaunay::Finite_facets_iterator& fit
) const
{
DynamicList<label> procsAttached(8);
const Cell_handle c1(fit->first);
const int oppositeVertex = fit->second;
const Cell_handle c2(c1->neighbor(oppositeVertex));
FixedList<label, 4> c1Procs(c1->processorsAttached());
FixedList<label, 4> c2Procs(c2->processorsAttached());
forAll(c1Procs, aPI)
{
if (findIndex(procsAttached, c1Procs[aPI] == -1))
{
procsAttached.append(c1Procs[aPI]);
}
if (findIndex(procsAttached, c2Procs[aPI] == -1))
{
procsAttached.append(c2Procs[aPI]);
}
}
return List<label>(procsAttached);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
#ifdef CGAL_INEXACT

View File

@ -423,10 +423,10 @@ void Foam::conformalVoronoiMesh::writeMesh
patches.setSize(nValidPatches);
// mesh.addFvPatches(patches);
mesh.addFvPatches(patches);
Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
mesh.addFvPatches(patches, false);
// Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
// mesh.addFvPatches(patches, false);
if (!mesh.write())
{
@ -456,7 +456,7 @@ void Foam::conformalVoronoiMesh::writeMesh
// cellCs.write();
// writeCellSizes(mesh);
writeCellSizes(mesh);
findRemainingProtrusionSet(mesh);
}

View File

@ -36,7 +36,10 @@ Description
#include <CGAL/Triangulation_3.h>
#include <CGAL/Triangulation_cell_base_with_circumcenter_3.h>
#include "indexedVertex.H"
#include "List.H"
#include "globalIndex.H"
#include "Pstream.H"
#include "Swap.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -187,7 +190,7 @@ public:
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const
{
if (!parallelDualVertex)
if (!parallelDualVertex())
{
return -1;
}
@ -212,6 +215,82 @@ public:
}
inline Foam::FixedList<Foam::label, 4> processorsAttached() const
{
if (!parallelDualVertex())
{
return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
}
Foam::FixedList<Foam::label, 4> procsAttached
(
Foam::Pstream::myProcNo()
);
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
procsAttached[i] = this->vertex(i)->procIndex();
}
}
return procsAttached;
}
// Using the globalIndex object, return a list of four (sorted) global
// Delaunay vertex indices that uniquely identify this tet in parallel
inline Foam::FixedList<Foam::label, 4> vertexGlobalIndices
(
const Foam::globalIndex& globalDelaunayVertexIndices
) const
{
// tetVertexGlobalIndices
Foam::FixedList<Foam::label, 4> tVGI(-1);
for (int i = 0; i < 4; i++)
{
Vertex_handle v = this->vertex(i);
// Finding the global index of each Delaunay vertex
if (v->referred())
{
// Referred vertices may have negative indices
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
v->procIndex(),
Foam::mag(v->index())
);
}
else
{
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
Foam::Pstream::myProcNo(),
v->index()
);
}
}
// bubble sort
for (int i = 0; i < tVGI.size(); i++)
{
for (int j = tVGI.size() - 1 ; j > i; j--)
{
if (tVGI[j - 1] > tVGI[j])
{
Foam::Swap(tVGI[j - 1], tVGI[j]);
}
}
}
return tVGI;
}
// Is the Delaunay cell part of the final dual mesh, i.e. any vertex form
// part of the internal or boundary definition
inline bool internalOrBoundaryDualVertex() const