Commit da64ab17 authored by franjo_j@hotmail.com's avatar franjo_j@hotmail.com
Browse files

Finished templated frontal marking. Usage for checking independent

regions and for removing domains.


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@29 fdcce57e-7e00-11e2-b579-49867b4cea03
parent d4617002
......@@ -36,8 +36,6 @@ Description
#include <omp.h>
# endif
#define DEBUGFrontalMarking
namespace Foam
{
......@@ -96,6 +94,61 @@ void frontalMarking
}
}
class graphNeiOp
{
// Private data
//- const reference to VRWGraph
const VRWGraph& neiGroups_;
public:
// Constructors
//- Construct from VRWGraph
inline graphNeiOp(const VRWGraph& neiGroups)
:
neiGroups_(neiGroups)
{}
// Public member functions
//- return the size of the graph
inline label size() const
{
return neiGroups_.size();
}
//- operator for finding neighbours
inline void operator()(const label groupI, DynList<label>& ng) const
{
ng = neiGroups_[groupI];
}
};
class graphSelectorOp
{
// Private data
//- const reference to VRWGraph
const VRWGraph& neiGroups_;
public:
// Constructors
//- Construct from VRWGraph
inline graphSelectorOp(const VRWGraph& neiGroups)
:
neiGroups_(neiGroups)
{}
// Public member functions
//- operator for selecting elements
inline bool operator()(const label groupI) const
{
if( (groupI < 0) || (groupI >= neiGroups_.size()) )
return false;
return true;
}
};
template<class labelListType, class neiOp, class filterOp>
label groupMarking
(
......@@ -226,89 +279,177 @@ label groupMarking
}
}
//- start processing connections between the group and merge the connected
//- ones into a new group
DynList<label> globalGroupLabel;
globalGroupLabel.setSize(nGroups);
globalGroupLabel = -1;
bool changed;
//- reduce the information about the groups
label counter(0);
do
forAll(neighbouringGroups, groupI)
{
changed = false;
if( globalGroupLabel[groupI] != -1 )
continue;
DynList<label> newGroupLabel;
newGroupLabel.setSize(globalGroupLabel.size());
newGroupLabel = -1;
DynList<label> connectedGroups;
frontalMarking
(
connectedGroups,
groupI,
graphNeiOp(neighbouringGroups),
graphSelectorOp(neighbouringGroups)
);
//- reduce the information about the groups
forAllReverse(neighbouringGroups, groupI)
{
forAllRow(neighbouringGroups, groupI, ngI)
{
const label neiGroup = neighbouringGroups(groupI, ngI);
forAll(connectedGroups, gI)
globalGroupLabel[connectedGroups[gI]] = counter;
if( neiGroup >= groupI )
continue;
++counter;
}
forAllRow(neighbouringGroups, groupI, i)
{
const label grI = neighbouringGroups(groupI, i);
neighbouringGroups.appendIfNotIn(neiGroup, grI);
}
}
}
nGroups = counter;
label counter(0);
forAll(neighbouringGroups, groupI)
{
if( newGroupLabel[groupI] != -1 )
continue;
forAll(neighbouringGroups, groupI)
{
if( globalGroupLabel[groupI] != -1 )
continue;
forAllRow(neighbouringGroups, groupI, ngI)
globalGroupLabel[neighbouringGroups(groupI, ngI)] = counter;
++counter;
}
if( Pstream::parRun() )
{
//- reduce the groups over processors of an MPI run
//- count the total number of groups over all processors
labelList nGroupsAtProc(Pstream::nProcs());
nGroupsAtProc[Pstream::myProcNo()] = nGroups;
Pstream::gatherList(nGroupsAtProc);
Pstream::scatterList(nGroupsAtProc);
forAllRow(neighbouringGroups, groupI, ngI)
newGroupLabel[neighbouringGroups(groupI, ngI)] = counter;
label startGroup(0), totalNumGroups(0);
for(label procI=0;procI<Pstream::nProcs();++procI)
{
totalNumGroups += nGroupsAtProc[procI];
++counter;
if( procI < Pstream::myProcNo() )
startGroup += nGroupsAtProc[procI];
}
if( counter < nGroups )
//- translate group labels
forAll(globalGroupLabel, groupI)
globalGroupLabel[groupI] += startGroup;
//- find the neighbouring groups
//- collect groups on other processors
//- this operator implements the algorithm for exchanging data
//- over processors and collects information which groups
//- are connected over inter-processor boundaries
std::map<label, DynList<label> > neiGroups;
neighbourCalculator.collectGroups
(
neiGroups,
elementInGroup,
globalGroupLabel
);
//- create a graph of connections
List<List<labelPair> > globalNeiGroups(Pstream::nProcs());
DynList<labelPair> connsAtProc;
for
(
std::map<label, DynList<label> >::const_iterator it =
neiGroups.begin();
it!=neiGroups.end();
++it
)
{
changed = true;
nGroups = counter;
const DynList<label>& ng = it->second;
forAll(ng, i)
connsAtProc.append(labelPair(it->first, ng[i]));
}
globalGroupLabel = newGroupLabel;
//- copy the connections into the global neighbour list
globalNeiGroups[Pstream::myProcNo()].setSize(connsAtProc.size());
forAll(connsAtProc, i)
globalNeiGroups[Pstream::myProcNo()][i] = connsAtProc[i];
//- communicate partial graphs to the master processor
Pstream::gatherList(globalNeiGroups);
if( Pstream::parRun() )
labelList allGroupsNewLabel;
if( Pstream::master() )
{
//- reduce the groups over processors
labelList nGroupsAtProc(Pstream::nProcs());
nGroupsAtProc[Pstream::myProcNo()] = nGroups;
//- collect the graph of connections for the whole system
VRWGraph allGroups(totalNumGroups);
forAll(allGroups, i)
allGroups[i].append(i);
Pstream::gatherList(nGroupsAtProc);
Pstream::scatterList(nGroupsAtProc);
forAll(globalNeiGroups, procI)
{
const List<labelPair>& connections = globalNeiGroups[procI];
label startGroup(0);
for(label procI=0;procI<Pstream::myProcNo();++procI)
startGroup += nGroupsAtProc[procI];
forAll(connections, i)
{
const labelPair& lp = connections[i];
//- find the neighbouring groups
std::map<label, DynList<label> > neiGroups;
allGroups.appendIfNotIn(lp.first(), lp.second());
allGroups.appendIfNotIn(lp.second(), lp.first());
}
}
//- collect groups on other processors
//- this operator implements the algorithm for exchanging data
//- over processors and collects information which groups
//- are connected ovr inter-processor boundaries
neighbourCalculator.collectGroups
(
neiGroups,
elementInGroup,
newGroupLabel,
startGroup
);
//- assign a global label to each group
allGroupsNewLabel.setSize(totalNumGroups);
allGroupsNewLabel = -1;
counter = 0;
forAll(allGroups, groupI)
{
if( allGroupsNewLabel[groupI] != -1 )
continue;
DynList<label> connectedGroups;
frontalMarking
(
connectedGroups,
groupI,
graphNeiOp(allGroups),
graphSelectorOp(allGroups)
);
forAll(connectedGroups, gI)
allGroupsNewLabel[connectedGroups[gI]] = counter;
++counter;
}
nGroups = counter;
}
} while( changed );
//- broadcast group labels from the master to other processors
Pstream::scatter(nGroups);
Pstream::scatter(allGroupsNewLabel);
//- assign correct group labels
forAll(globalGroupLabel, groupI)
{
const label ngI = globalGroupLabel[groupI];
globalGroupLabel[groupI] = allGroupsNewLabel[ngI];
}
}
//- set the global group label
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 50)
# endif
forAll(elementInGroup, elI)
{
if( elementInGroup[elI] < 0 )
......
......@@ -38,6 +38,8 @@ Description
#include <omp.h>
# endif
//#define DEBUGRemoveDomains
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -45,6 +47,141 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class meshNeighbourOperator
{
const polyMeshGen& mesh_;
public:
meshNeighbourOperator(const polyMeshGen& mesh)
:
mesh_(mesh)
{}
label size() const
{
return mesh_.cells().size();
}
void operator()(const label cellI, DynList<label>& neighbourCells) const
{
neighbourCells.clear();
const labelList& owner = mesh_.owner();
const labelList& neighbour = mesh_.neighbour();
const cell& c = mesh_.cells()[cellI];
forAll(c, fI)
{
label nei = owner[c[fI]];
if( nei == cellI )
nei = neighbour[c[fI]];
if( nei >= 0 )
neighbourCells.append(nei);
}
}
template<class labelListType>
void collectGroups
(
std::map<label, DynList<label> >& neiGroups,
const labelListType& elementInGroup,
const DynList<label>& localGroupLabel
) const
{
const PtrList<writeProcessorPatch>& procBoundaries =
mesh_.procBoundaries();
const labelList& owner = mesh_.owner();
//- send the data to other processors
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label size = procBoundaries[patchI].patchSize();
labelList groupOwner(procBoundaries[patchI].patchSize());
for(label faceI=0;faceI<size;++faceI)
{
const label groupI = elementInGroup[owner[start+faceI]];
if( groupI < 0 )
{
groupOwner[faceI] = -1;
continue;
}
groupOwner[faceI] = localGroupLabel[groupI];
}
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
groupOwner.byteSize()
);
toOtherProc << groupOwner;
}
//- receive data from other processors
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
labelList receivedData;
IPstream fromOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
fromOtherProc >> receivedData;
forAll(receivedData, faceI)
{
if( receivedData[faceI] < 0 )
continue;
const label groupI = elementInGroup[owner[start+faceI]];
if( groupI < 0 )
continue;
DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
//- store the connection over the inter-processor boundary
ng.appendIfNotIn(receivedData[faceI]);
}
}
}
};
class meshSelectorOperator
{
const LongList<DynList<label, 2> >& intersectedCells_;
public:
meshSelectorOperator(const LongList<DynList<label, 2> >& intersectedCells)
:
intersectedCells_(intersectedCells)
{}
bool operator()(const label cellI) const
{
if( intersectedCells_[cellI].size() )
return false;
return true;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
const VRWGraph& removeCellsInSelectedDomains::cellsIntersectedBySurfaceFacets()
{
if( !surfIntersectionPtr_ )
......@@ -172,71 +309,22 @@ void removeCellsInSelectedDomains::findAndRemoveCells()
{
const DynList<label, 2>& doms = intersectedCells[i];
forAll(doms, j)
mesh_.addCellToSubset(domainIds[doms[j]], i);
forAll(doms, ii)
mesh_.addCellToSubset(domainIds[doms[i]], i);
};
# endif
//- TODO: implement the group marking algorithm properly using templates
//- find islands of cells which are not intersected by the selected domains
const cellListPMG& cells = mesh_.cells();
const labelList& owner = mesh_.owner();
const labelList& neighbour = mesh_.neighbour();
labelListPMG findCellGroups(intersectedCells.size(), -1);
label nGroups(0);
forAll(findCellGroups, cellI)
{
if( findCellGroups[cellI] != -1 )
continue;
if( intersectedCells[cellI].size() )
continue;
findCellGroups[cellI] = nGroups;
//const cellListPMG& cells = mesh_.cells();
labelListPMG findCellGroups;
labelListPMG front;
front.append(cellI);
while( front.size() )
{
const label cLabel = front.removeLastElement();
const cell& c = cells[cLabel];
//- find neighbouring elements
DynList<label> neis;
forAll(c, fI)
{
const label faceI = c[fI];
if( owner[faceI] == cLabel && neighbour[faceI] >= 0 )
{
neis.append(neighbour[faceI]);
}
else if( neighbour[faceI] == cLabel )
{
neis.append(owner[faceI]);
}
}
//- check which neighbours pass the criteria
forAll(neis, i)
{
const label nei = neis[i];
if( findCellGroups[nei] != -1 )
continue;
if( intersectedCells[nei].size() )
continue;
findCellGroups[nei] = nGroups;
front.append(nei);
}
}
++nGroups;
}
const label nGroups =
help::groupMarking
(
findCellGroups,
meshNeighbourOperator(mesh_),
meshSelectorOperator(intersectedCells)
);
# ifdef DEBUGRemoveDomains
Info << "Number of groups " << nGroups << endl;
......@@ -254,6 +342,9 @@ void removeCellsInSelectedDomains::findAndRemoveCells()
//- collect which domains are assigned to facets neighbouring
//- groups
const labelList& owner = mesh_.owner();
const labelList& neighbour = mesh_.neighbour();
List<DynList<label> > neiDomains(nGroups);
for(label faceI=0;faceI<mesh_.nInternalFaces();++faceI)
{
......@@ -278,6 +369,60 @@ void removeCellsInSelectedDomains::findAndRemoveCells()
}
}
if( Pstream::parRun() )
{
const PtrList<writeProcessorPatch>& procBoundaries =
mesh_.procBoundaries();
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label size = procBoundaries[patchI].patchSize();
labelList sendGroup(size);
for(label fI=0;fI<size;++fI)
sendGroup[fI] = findCellGroups[owner[start+fI]];
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
sendGroup.byteSize()
);
toOtherProc << sendGroup;
}
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
labelList receivedGroup;