/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | cfMesh: A library for mesh generation
\\ / O peration |
\\ / A nd | Author: Franjo Juretic (franjo.juretic@c-fields.com)
\\/ M anipulation | Copyright (C) Creative Fields, Ltd.
-------------------------------------------------------------------------------
License
This file is part of cfMesh.
cfMesh 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.
cfMesh 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 cfMesh. If not, see .
Description
\*---------------------------------------------------------------------------*/
#include "meshSurfaceEngine.H"
#include "demandDrivenData.H"
#include "boolList.H"
#include "helperFunctions.H"
#include "VRWGraphSMPModifier.H"
#include "labelledPoint.H"
#include "HashSet.H"
#include
#include
# ifdef USE_OMP
#include
# endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void meshSurfaceEngine::calculateBoundaryFaces() const
{
if( mesh_.boundaries().size() != 0 )
{
const faceListPMG& faces = mesh_.faces();
const PtrList& boundaries = mesh_.boundaries();
label nBoundaryFaces(0);
if( activePatch_ < 0 )
{
//- take all patches
forAll(boundaries, patchI)
nBoundaryFaces += boundaries[patchI].patchSize();
boundaryFacesPtr_ =
new faceList::subList
(
faces,
nBoundaryFaces,
boundaries[0].patchStart()
);
}
else if( activePatch_ < boundaries.size() )
{
nBoundaryFaces = boundaries[activePatch_].patchSize();
boundaryFacesPtr_ =
new faceList::subList
(
faces,
nBoundaryFaces,
boundaries[activePatch_].patchStart()
);
}
else
{
FatalErrorIn
(
"void meshSurfaceEngine::calculateBoundaryFaces() const"
) << "Cannot select boundary faces. Invalid patch index "
<< activePatch_ << exit(FatalError);
}
reduce(nBoundaryFaces, sumOp());
Info << "Found " << nBoundaryFaces << " boundary faces " << endl;
}
else
{
FatalErrorIn
(
"void meshSurfaceEngine::calculateBoundaryFaces() const"
) << "Boundary faces are not at the end of the face list!"
<< exit(FatalError);
}
}
void meshSurfaceEngine::calculateBoundaryOwners() const
{
const labelList& owner = mesh_.owner();
const faceList::subList& boundaryFaces = this->boundaryFaces();
if( !boundaryFaceOwnersPtr_ )
boundaryFaceOwnersPtr_ = new labelList(boundaryFaces.size());
labelList& owners = *boundaryFaceOwnersPtr_;
const label start = mesh_.boundaries()[0].patchStart();
# ifdef USE_OMP
# pragma omp parallel for schedule(static, 1)
# endif
forAll(boundaryFaces, fI)
owners[fI] = owner[start+fI];
}
void meshSurfaceEngine::calculateBoundaryNodes() const
{
//- mark boundary points
label pointI(0);
if( !bppPtr_ )
bppPtr_ = new labelList(mesh_.points().size(), -1);
labelList& bp = *bppPtr_;
const faceList::subList& boundaryFaces = this->boundaryFaces();
boolList isBndPoint(bp.size(), false);
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# pragma omp parallel for num_threads(nThreads) schedule(static, 1)
# endif
forAll(boundaryFaces, bfI)
{
const face& bf = boundaryFaces[bfI];
forAll(bf, pI)
isBndPoint[bf[pI]] = true;
}
forAll(isBndPoint, pI)
{
if( isBndPoint[pI] )
bp[pI] = pointI++;
}
if( Pstream::parRun() )
{
const faceListPMG& faces = mesh_.faces();
const PtrList& procBoundaries =
mesh_.procBoundaries();
//- exchange information with processors
//- this is need sometimes to find all nodes at the boundary
bool found;
do
{
found = false;
//- send bnd nodes to other processor
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label end = start + procBoundaries[patchI].patchSize();
labelLongList dts;
labelHashSet addedPoint;
for(label faceI=start;faceI> receiveData;
label counter(0);
while( counter < receiveData.size() )
{
const label fI = receiveData[counter++];
const label pI = receiveData[counter++];
if( bp[faces[start+fI][pI]] == -1 )
{
bp[faces[start+fI][pI]] = pointI++;
found = true;
}
}
}
reduce(found, maxOp());
} while( found );
}
if( !boundaryPointsPtr_ )
boundaryPointsPtr_ = new labelList();
labelList& boundaryPoints = *boundaryPointsPtr_;
boundaryPoints.setSize(pointI);
//- fill the boundaryPoints list
# ifdef USE_OMP
# pragma omp parallel for num_threads(nThreads) schedule(static, 1)
# endif
forAll(bp, bpI)
{
if( bp[bpI] != -1 )
boundaryPoints[bp[bpI]] = bpI;
}
}
void meshSurfaceEngine::calculateBoundaryFacePatches() const
{
const faceList::subList& bFaces = this->boundaryFaces();
boundaryFacePatchPtr_ = new labelList(bFaces.size());
labelList& facePatch = *boundaryFacePatchPtr_;
label faceI(0);
const PtrList& boundaries = mesh_.boundaries();
forAll(boundaries, patchI)
{
const label nFaces = boundaries[patchI].patchSize();
for(label patchFaceI=0;patchFaceIboundaryPoints();
const faceList::subList& bFaces = this->boundaryFaces();
//- create boundary points
const labelList& bp = this->bp();
labelLongList npf;
# ifdef USE_OMP
label nThreads = 3 * omp_get_num_procs();
if( bPoints.size() < 1000 )
nThreads = 1;
# else
const label nThreads(1);
# endif
label minRow(INT_MAX), maxRow(0);
List > > dataForOtherThreads(nThreads);
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
const label threadI = omp_get_thread_num();
# else
const label threadI(0);
# endif
List >& dot = dataForOtherThreads[threadI];
dot.setSize(nThreads);
//- find min and max entry in the graph
//- they are used for assigning ranges of values local for each process
label localMinRow(minRow), localMaxRow(0);
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(bFaces, bfI)
{
const face& bf = bFaces[bfI];
forAll(bf, pI)
{
const label bpI = bp[bf[pI]];
localMaxRow = Foam::max(localMaxRow, bpI);
localMinRow = Foam::min(localMinRow, bpI);
}
}
++localMaxRow;
# ifdef USE_OMP
# pragma omp critical
# endif
{
minRow = Foam::min(minRow, localMinRow);
minRow = Foam::max(minRow, 0);
maxRow = Foam::max(maxRow, localMaxRow);
npf.setSize(maxRow);
}
# ifdef USE_OMP
# pragma omp barrier
# endif
//- initialise appearances
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
for(label i=0;i& data =
dataForOtherThreads[i][threadI];
forAll(data, j)
++npf[data[j].first()];
}
# ifdef USE_OMP
# pragma omp barrier
# endif
//- allocate graph
# ifdef USE_OMP
# pragma omp master
# endif
{
VRWGraphSMPModifier(pointFacesAddr).setSizeAndRowSize(npf);
VRWGraphSMPModifier(pointInFaceAddr).setSizeAndRowSize(npf);
}
# ifdef USE_OMP
# pragma omp barrier
# endif
for(label i=localMin;i& data =
dataForOtherThreads[i][threadI];
forAll(data, j)
{
const label bpI = data[j].first();
const label bfI = data[j].second();
pointFacesAddr(bpI, npf[bpI]) = bfI;
pointInFaceAddr(bpI, npf[bpI]) =
bFaces[bfI].which(bPoints[bpI]);
++npf[bpI];
}
}
//- update data local to the processor
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(bFaces, bfI)
{
const face& bf = bFaces[bfI];
forAll(bf, pI)
{
const label bpI = bp[bf[pI]];
if( (bpI >= localMin) && (bpI < localMax) )
{
pointInFaceAddr(bpI, npf[bpI]) = pI;
pointFacesAddr(bpI, npf[bpI]++) = bfI;
}
}
}
//- update data from the processors with higher labels
for(label i=threadI+1;i& data =
dataForOtherThreads[i][threadI];
forAll(data, j)
{
const label bpI = data[j].first();
const label bfI = data[j].second();
pointFacesAddr(bpI, npf[bpI]) = bfI;
pointInFaceAddr(bpI, npf[bpI]) =
bFaces[bfI].which(bPoints[bpI]);
++npf[bpI];
}
}
}
pointFacesAddr.setSize(bPoints.size());
pointInFaceAddr.setSize(bPoints.size());
}
void meshSurfaceEngine::calculatePointPatches() const
{
if( !pointPatchesPtr_ )
pointPatchesPtr_ = new VRWGraph();
VRWGraph& pPatches = *pointPatchesPtr_;
const labelList& facePatch = boundaryFacePatches();
const VRWGraph& pFaces = pointFaces();
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# endif
labelList npPatches(pFaces.size());
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(npPatches, i)
npPatches[i] = 0;
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(pFaces, bpI)
{
DynList pf;
forAllRow(pFaces, bpI, pfI)
pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);
npPatches[bpI] = pf.size();
}
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
# endif
VRWGraphSMPModifier(pPatches).setSizeAndRowSize(npPatches);
# ifdef USE_OMP
# pragma omp barrier
# pragma omp for schedule(static)
# endif
forAll(pFaces, bpI)
{
DynList pf;
forAllRow(pFaces, bpI, pfI)
pf.appendIfNotIn(facePatch[pFaces(bpI, pfI)]);
pPatches.setRow(bpI, pf);
}
}
if( Pstream::parRun() )
{
const labelList& globalPointLabel = globalBoundaryPointLabel();
const VRWGraph& bpAtProcs = this->bpAtProcs();
const Map& globalToLocal = globalToLocalBndPointAddressing();
std::map exchangeData;
forAllConstIter(Map, globalToLocal, iter)
{
const label bpI = iter();
forAllRow(bpAtProcs, bpI, procI)
{
const label neiProc = bpAtProcs(bpI, procI);
if( neiProc == Pstream::myProcNo() )
continue;
if( exchangeData.find(neiProc) == exchangeData.end() )
{
exchangeData.insert
(
std::make_pair(neiProc, labelLongList())
);
}
labelLongList& dataToSend = exchangeData[neiProc];
//- prepare data which will be sent
//- data is sent as follows
//- 1. global point label
//- 2. number of local patches for point
//- 3. patch labels for a given point
dataToSend.append(globalPointLabel[bpI]);
dataToSend.append(pPatches.sizeOfRow(bpI));
forAllRow(pPatches, bpI, patchI)
dataToSend.append(pPatches(bpI, patchI));
}
}
//- exchange data with other processors
labelLongList receivedData;
help::exchangeMap(exchangeData, receivedData);
label counter(0);
while( counter < receivedData.size() )
{
const label bpI = globalToLocal[receivedData[counter++]];
const label nPatches = receivedData[counter++];
for(label i=0;iboundaryPoints();
const faceList::subList& bFaces = this->boundaryFaces();
const VRWGraph& pFaces = this->pointFaces();
const labelList& bp = this->bp();
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# endif
labelList npp(boundaryPoints.size());
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(npp, i)
npp[i] = 0;
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(pFaces, bpI)
{
DynList pPoints;
forAllRow(pFaces, bpI, pfI)
{
const face& bf = bFaces[pFaces(bpI, pfI)];
const label pos = bf.which(boundaryPoints[bpI]);
pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
}
npp[bpI] = pPoints.size();
}
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
# endif
VRWGraphSMPModifier(pointPoints).setSizeAndRowSize(npp);
# ifdef USE_OMP
# pragma omp barrier
# pragma omp for schedule(static)
# endif
forAll(pFaces, bpI)
{
DynList pPoints;
forAllRow(pFaces, bpI, pfI)
{
const face& bf = bFaces[pFaces(bpI, pfI)];
const label pos = bf.which(boundaryPoints[bpI]);
pPoints.appendIfNotIn(bp[bf.nextLabel(pos)]);
pPoints.appendIfNotIn(bp[bf.prevLabel(pos)]);
}
pointPoints.setRow(bpI, pPoints);
}
}
if( Pstream::parRun() )
{
//- this is needed to make the connection matrix symmetric
//- on all processors. In some cases the points on a given processor
//- may not be connected because of a single layer of faces on some
//- other processor. P0, P0 | P1 | P0 P0
const labelList& globalPointLabel = this->globalBoundaryPointLabel();
const Map& globalToLocal =
this->globalToLocalBndPointAddressing();
const VRWGraph& bpAtProcs = this->bpAtProcs();
std::map exchangeData;
forAllConstIter(Map, globalToLocal, iter)
{
const label bpI = iter();
DynList neiToSend;
forAllRow(pointPoints, bpI, j)
{
const label bpJ = pointPoints(bpI, j);
if( bpAtProcs.sizeOfRow(bpJ) != 0 )
neiToSend.append(bpJ);
}
forAllRow(bpAtProcs, bpI, procI)
{
const label neiProc = bpAtProcs(bpI, procI);
if( neiProc == Pstream::myProcNo() )
continue;
if( exchangeData.find(neiProc) == exchangeData.end() )
exchangeData.insert(std::make_pair(neiProc,labelLongList()));
if( neiToSend.size() != 0 )
{
labelLongList& dts = exchangeData[neiProc];
dts.append(globalPointLabel[bpI]);
dts.append(neiToSend.size());
forAll(neiToSend, i)
dts.append(globalPointLabel[neiToSend[i]]);
}
}
}
labelLongList receivedData;
help::exchangeMap(exchangeData, receivedData);
label counter(0);
while( counter < receivedData.size() )
{
const label bpI = globalToLocal[receivedData[counter++]];
const label size = receivedData[counter++];
for(label i=0;i 1000 ) schedule(dynamic, 50)
# endif
forAll(pFaces, pI)
{
vector normal(vector::zero);
forAllRow(pFaces, pI, pfI)
normal += fNormals[pFaces(pI, pfI)];
const scalar d = mag(normal);
if( d > VSMALL )
{
normal /= d;
}
else
{
normal = vector::zero;
}
(*pointNormalsPtr_)[pI] = normal;
}
updatePointNormalsAtProcBoundaries();
}
void meshSurfaceEngine::calculateFaceNormals() const
{
const faceList::subList& bFaces = this->boundaryFaces();
const pointFieldPMG& points = mesh_.points();
faceNormalsPtr_ = new vectorField(bFaces.size());
# ifdef USE_OMP
# pragma omp parallel for if( bFaces.size() > 1000 )
# endif
forAll(bFaces, bfI)
{
const face& bf = bFaces[bfI];
faceNormalsPtr_->operator[](bfI) = bf.normal(points);
}
}
void meshSurfaceEngine::calculateFaceCentres() const
{
const faceList::subList& bFaces = this->boundaryFaces();
const pointFieldPMG& points = mesh_.points();
faceCentresPtr_ = new vectorField(bFaces.size());
# ifdef USE_OMP
# pragma omp parallel for if( bFaces.size() > 1000 )
# endif
forAll(bFaces, bfI)
faceCentresPtr_->operator[](bfI) = bFaces[bfI].centre(points);
}
void meshSurfaceEngine::updatePointNormalsAtProcBoundaries() const
{
if( !Pstream::parRun() )
return;
const VRWGraph& pFaces = pointFaces();
const vectorField& fNormals = faceNormals();
const labelList& globalPointLabel = this->globalBoundaryPointLabel();
const Map& globalToLocal =
this->globalToLocalBndPointAddressing();
const VRWGraph& bpAtProcs = this->bpAtProcs();
vectorField& pNormals = *pointNormalsPtr_;
//- create data which will be sent to other processors
std::map > exchangeData;
forAllConstIter(Map, globalToLocal, iter)
{
const label bpI = iter();
vector& n = pNormals[bpI];
n = vector::zero;
forAllRow(pFaces, bpI, pfI)
n += fNormals[pFaces(bpI, pfI)];
forAllRow(bpAtProcs, bpI, procI)
{
const label neiProc = bpAtProcs(bpI, procI);
if( neiProc == Pstream::myProcNo() )
continue;
if( exchangeData.find(neiProc) == exchangeData.end() )
{
exchangeData.insert
(
std::make_pair(neiProc, LongList())
);
}
exchangeData[neiProc].append
(
labelledPoint(globalPointLabel[bpI], n)
);
}
}
//- exchange data with other procs
LongList receivedData;
help::exchangeMap(exchangeData, receivedData);
forAll(receivedData, i)
{
const label bpI = globalToLocal[receivedData[i].pointLabel()];
pNormals[bpI] += receivedData[i].coordinates();
}
//- normalize vectors
# ifdef USE_OMP
# pragma omp parallel for if( bpAtProcs.size() > 1000 ) \
schedule(guided)
# endif
forAll(bpAtProcs, bpI)
{
if( bpAtProcs.sizeOfRow(bpI) == 0 )
continue;
vector normal = pNormals[bpI];
const scalar d = mag(normal);
if( d > VSMALL )
{
normal /= d;
}
else
{
normal = vector::zero;
}
pNormals[bpI] = normal;
}
}
void meshSurfaceEngine::calculateEdgesAndAddressing() const
{
const VRWGraph& pFaces = pointFaces();
const faceList::subList& bFaces = boundaryFaces();
const labelList& bPoints = boundaryPoints();
const labelList& bp = this->bp();
edgesPtr_ = new edgeList();
edgeList& edges = *edgesPtr_;
bpEdgesPtr_ = new VRWGraph();
VRWGraph& bpEdges = *bpEdgesPtr_;
# ifdef USE_OMP
label nThreads = 3 * omp_get_num_procs();
if( pFaces.size() < 1000 )
nThreads = 1;
# else
const label nThreads(1);
# endif
labelList nEdgesForThread(nThreads);
label edgeI(0);
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
{
edgeLongList edgesHelper;
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(pFaces, bpI)
{
std::set > edgesAtPoint;
forAllRow(pFaces, bpI, pfI)
{
const label bfI = pFaces(bpI, pfI);
const face& bf = bFaces[bfI];
const label pos = bf.which(bPoints[bpI]);
if( bp[bf.nextLabel(pos)] >= bpI )
{
edgesAtPoint.insert
(
std::make_pair(bf[pos], bf.nextLabel(pos))
);
}
if( bp[bf.prevLabel(pos)] >= bpI )
{
edgesAtPoint.insert
(
std::make_pair(bf[pos], bf.prevLabel(pos))
);
}
}
std::set >::const_iterator it;
for(it=edgesAtPoint.begin();it!=edgesAtPoint.end();++it)
edgesHelper.append(edge(it->first, it->second));
}
//- this enables other threads to see the number of edges
//- generated by each thread
# ifdef USE_OMP
const label threadI = omp_get_thread_num();
# else
const label threadI(0);
# endif
nEdgesForThread[threadI] = edgesHelper.size();
# ifdef USE_OMP
# pragma omp critical
# endif
edgeI += edgesHelper.size();
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
# endif
edgesPtr_->setSize(edgeI);
# ifdef USE_OMP
# pragma omp barrier
# endif
//- find the starting position of the edges generated by this thread
//- in the global list of edges
label localStart(0);
for(label i=0;ioperator[](localStart+i) = edgesHelper[i];
}
//- set the bpEdges
VRWGraphSMPModifier(bpEdges).reverseAddressing(bp, edges);
bpEdges.setSize(pFaces.size());
if( !Pstream::parRun() )
return;
bool addEdges;
do
{
addEdges = false;
//- mark boundary edges for processors which do not contain
//- boundary faces. This procedure is needed to identify boundary
//- edges which are not part of any boundary face on their processor
const faceListPMG& faces = mesh_.faces();
const PtrList& procBoundaries =
mesh_.procBoundaries();
//- send boundary edges to neighbour processors
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label end = start + procBoundaries[patchI].patchSize();
labelLongList dts;
for(label faceI=start;faceI> receivedEdges;
const label start = procBoundaries[patchI].patchStart();
label nReceivedEdges(0);
while( nReceivedEdges < receivedEdges.size() )
{
const face& f = faces[start+receivedEdges[nReceivedEdges++]];
const label eI = receivedEdges[nReceivedEdges++];
const edge e = f.faceEdge(eI);
const label s = bp[e.start()];
bool found(false);
forAllRow(bpEdges, s, peI)
if( edges[bpEdges(s, peI)] == e )
{
found = true;
break;
}
if( !found )
{
//- create a new edge
addEdges = true;
edges.newElmt(edgeI) = e;
bpEdges.append(bp[e.start()], edgeI);
bpEdges.append(bp[e.end()], edgeI);
++edgeI;
}
}
}
reduce(addEdges, maxOp());
} while( addEdges );
edges.setSize(edgeI);
}
void meshSurfaceEngine::calculateFaceEdgesAddressing() const
{
const faceList::subList& bFaces = this->boundaryFaces();
const labelList& bp = this->bp();
const edgeList& edges = this->edges();
const VRWGraph& bpEdges = this->boundaryPointEdges();
faceEdgesPtr_ = new VRWGraph(bFaces.size());
VRWGraph& faceEdges = *faceEdgesPtr_;
labelList nfe(bFaces.size());
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(bFaces, bfI)
nfe[bfI] = bFaces[bfI].size();
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
{
# endif
VRWGraphSMPModifier(faceEdges).setSizeAndRowSize(nfe);
# ifdef USE_OMP
}
# pragma omp barrier
# pragma omp for schedule(dynamic, 100)
# endif
forAll(faceEdges, bfI)
{
const face& bf = bFaces[bfI];
forAll(bf, eI)
{
const edge e = bf.faceEdge(eI);
const label bps = bp[e.start()];
forAllRow(bpEdges, bps, peI)
{
const label beI = bpEdges(bps, peI);
const edge& ee = edges[beI];
if( e == ee )
{
faceEdges(bfI, eI) = beI;
break;
}
}
}
}
}
}
void meshSurfaceEngine::calculateEdgeFacesAddressing() const
{
const faceList::subList& bFaces = this->boundaryFaces();
const VRWGraph& pointFaces = this->pointFaces();
const edgeList& edges = this->edges();
const labelList& bp = this->bp();
edgeFacesPtr_ = new VRWGraph();
VRWGraph& edgeFaces = *edgeFacesPtr_;
labelList nef(edges.size());
# ifdef USE_OMP
const label nThreads = 3 * omp_get_num_procs();
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(nef, edgeI)
nef[edgeI] = 0;
# ifdef USE_OMP
# pragma omp for schedule(static)
# endif
forAll(edges, edgeI)
{
const edge& ee = edges[edgeI];
const label bpI = bp[ee.start()];
forAllRow(pointFaces, bpI, pfI)
{
const label bfI = pointFaces(bpI, pfI);
const face& bf = bFaces[bfI];
forAll(bf, eI)
{
if( bf.faceEdge(eI) == ee )
{
++nef[edgeI];
break;
}
}
}
}
# ifdef USE_OMP
# pragma omp barrier
# pragma omp master
# endif
VRWGraphSMPModifier(edgeFaces).setSizeAndRowSize(nef);
# ifdef USE_OMP
# pragma omp barrier
# pragma omp for schedule(static)
# endif
forAll(edges, edgeI)
{
const edge& ee = edges[edgeI];
const label bpI = bp[ee.start()];
//- find boundary faces attached to this edge
DynList eFaces;
forAllRow(pointFaces, bpI, pfI)
{
const label bfI = pointFaces(bpI, pfI);
const face& bf = bFaces[bfI];
forAll(bf, eI)
{
if( bf.faceEdge(eI) == ee )
{
eFaces.append(bfI);
break;
}
}
}
//- the face that owns the edge shall be the first one in the list
// TODO: find out whether this will be necessary
if( eFaces.size() == 2 )
{
const face& bf = bFaces[eFaces[1]];
const label pos = bf.which(ee.start());
if( bf.nextLabel(pos) == ee.end() )
{
//- this face shall be the first one in the list
const label helper = eFaces[0];
eFaces[0] = eFaces[1];
eFaces[1] = helper;
}
}
edgeFaces.setRow(edgeI, eFaces);
}
}
}
void meshSurfaceEngine::calculateEdgePatchesAddressing() const
{
edgePatchesPtr_ = new VRWGraph();
VRWGraph& edgePatches = *edgePatchesPtr_;
const VRWGraph& edgeFaces = this->edgeFaces();
const labelList& facePatch = this->boundaryFacePatches();
edgePatches.setSize(edgeFaces.size());
forAll(edgeFaces, eI)
{
DynList ePatches;
forAllRow(edgeFaces, eI, i)
{
const label patchI = facePatch[edgeFaces(eI, i)];
ePatches.appendIfNotIn(patchI);
}
edgePatches.setRow(eI, ePatches);
}
if( Pstream::parRun() )
{
const Map& globalToLocal = globalToLocalBndEdgeAddressing();
const Map& otherPatch = this->otherEdgeFacePatch();
forAllConstIter(Map, globalToLocal, it)
{
const label beI = it();
edgePatches.appendIfNotIn(beI, otherPatch[beI]);
}
}
}
void meshSurfaceEngine::calculateFaceFacesAddressing() const
{
const VRWGraph& edgeFaces = this->edgeFaces();
const faceList::subList& bFaces = boundaryFaces();
faceFacesPtr_ = new VRWGraph(bFaces.size());
VRWGraph& faceFaces = *faceFacesPtr_;
forAll(bFaces, bfI)
faceFaces.setRowSize(bfI, bFaces[bfI].size());
labelList nAppearances(bFaces.size(), 0);
forAll(edgeFaces, efI)
{
if( edgeFaces.sizeOfRow(efI) == 2 )
{
const label f0 = edgeFaces(efI, 0);
const label f1 = edgeFaces(efI, 1);
faceFaces(f0, nAppearances[f0]++) = f1;
faceFaces(f1, nAppearances[f1]++) = f0;
}
else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) == 1) )
{
const label f0 = edgeFaces(efI, 0);
faceFaces(f0, nAppearances[f0]++) = -1;
}
else if( Pstream::parRun() && (edgeFaces.sizeOfRow(efI) != 0 ) )
{
FatalErrorIn
(
"void meshSurfaceEngine::calculateFaceFacesAddressing() const"
) << "The surface of the mesh is invalid!"
<< " The number of faces containing edge " << efI
<< " is " << edgeFaces.sizeOfRow(efI)
<< " Cannot continue" << exit(FatalError);
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //