partTriMeshParallelAddressing.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 <http://www.gnu.org/licenses/>.
 
Description
 
\*---------------------------------------------------------------------------*/
 
#include "demandDrivenData.H"
#include "meshSurfacePartitioner.H"
#include "partTriMesh.H"
#include "triSurfModifier.H"
#include "meshSurfaceEngine.H"
#include "helperFunctionsPar.H"
#include "parTriFace.H"
 
#include <map>
 
//#define DEBUGSmooth
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
void partTriMesh::createParallelAddressing
(
    const labelList& nodeLabelForPoint,
    const labelList& /*nodeLabelForFace*/
)
{
    const meshSurfaceEngine& mse = mPart_.surfaceEngine();
 
    const pointField& pts = surf_.points();
 
    //- vertices marked as SMOOTH are used by the smoother
    const direction useType = SMOOTH;
 
    //- allocate global point labels
    if( !globalPointLabelPtr_ )
        globalPointLabelPtr_ = new labelLongList();
    labelLongList& globalPointLabel = *globalPointLabelPtr_;
    globalPointLabel.setSize(pts.size());
    globalPointLabel = -1;
 
    //- allocated point-processors addressing
    if( !pAtProcsPtr_ )
        pAtProcsPtr_ = new VRWGraph();
    VRWGraph& pProcs = *pAtProcsPtr_;
    pProcs.setSize(0);
    pProcs.setSize(pts.size());
 
    //- allocate global-to-local point addressing
    if( !globalToLocalPointAddressingPtr_ )
        globalToLocalPointAddressingPtr_ = new Map<label>();
    Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
    globalToLocal.clear();
 
    //- allocate storage for points at parallel boundaries
    if( !pAtParallelBoundariesPtr_ )
        pAtParallelBoundariesPtr_ = new labelLongList();
    labelLongList& pAtParallelBoundaries = *pAtParallelBoundariesPtr_;
    pAtParallelBoundaries.clear();
 
    //- create point-processors addressing
    std::map<label, labelLongList> exchangeData;
    std::map<label, labelLongList>::iterator iter;
 
    const Map<label>& globalToLocalPointAddressing =
        mse.globalToLocalBndPointAddressing();
    const VRWGraph& pAtProcs = mse.bpAtProcs();
    const DynList<label>& pNeiProcs = mse.bpNeiProcs();
 
    forAll(pNeiProcs, procI)
        exchangeData.insert(std::make_pair(pNeiProcs[procI], labelLongList()));
 
    //- make sure that the same vertices are marked for smoothing on all procs
    //- this is performed by sending the labels of vertices which are not used
    //- for tet mesh creation and the tet mesh vertices which are not moved
    forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
    {
        const label pI = it();
 
        if(
            nodeLabelForPoint[pI] == -1 ||
            !pointType_[nodeLabelForPoint[pI]]
        )
        {
            forAllRow(pAtProcs, pI, procI)
            {
                const label neiProc = pAtProcs(pI, procI);
                if( neiProc == Pstream::myProcNo() )
                    continue;
 
                exchangeData[neiProc].append(it.key());
            }
        }
    }
 
    //- exchange data with other processors
    labelLongList receivedData;
    help::exchangeMap(exchangeData, receivedData);
 
    //- set the values according to other processors
    forAll(receivedData, i)
    {
        const label pointI = globalToLocalPointAddressing[receivedData[i]];
 
        if( nodeLabelForPoint[pointI] == -1 )
            continue;
 
        pointType_[nodeLabelForPoint[pointI]] = NONE;
    }
 
    for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
        iter->second.clear();
 
    //- start creating global-to-local addressing
    //- find the starting point labels
    label startPoint(0), nLocalPoints(0), nSharedPoints(0);
 
    //- count the number of points at processor boundaries
    forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
    {
        const label pI = it();
 
        if( nodeLabelForPoint[pI] == -1 )
            continue;
        if( !(pointType_[nodeLabelForPoint[pI]] & useType) )
            continue;
 
        ++nSharedPoints;
 
        label pMin(Pstream::myProcNo());
        forAllRow(pAtProcs, pI, procI)
            pMin = Foam::min(pMin, pAtProcs(pI, procI));
 
        if( pMin == Pstream::myProcNo() )
            ++nLocalPoints;
    }
 
    labelList nPointsAtProc(Pstream::nProcs());
    nSharedPoints -= nLocalPoints;
    nPointsAtProc[Pstream::myProcNo()] = pts.size() - nSharedPoints;
    Pstream::gatherList(nPointsAtProc);
    Pstream::scatterList(nPointsAtProc);
 
    for(label i=0;i<Pstream::myProcNo();++i)
        startPoint += nPointsAtProc[i];
 
    //- create global labels for points at processor boundaries
    forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
    {
        const label pI = it();
 
        if( nodeLabelForPoint[pI] == -1 )
            continue;
 
        const label pLabel = nodeLabelForPoint[pI];
 
        if( !(pointType_[pLabel] & useType) )
            continue;
 
        label pMin(Pstream::myProcNo());
        forAllRow(pAtProcs, pI, procI)
        {
            const label neiProc = pAtProcs(pI, procI);
            pProcs.append(pLabel, neiProc);
            pMin = Foam::min(pMin, neiProc);
        }
 
        if( pMin != Pstream::myProcNo() )
            continue;
 
        globalPointLabel[pLabel] = startPoint++;
 
        forAllRow(pAtProcs, pI, procI)
        {
            const label neiProc = pAtProcs(pI, procI);
 
            if( neiProc == Pstream::myProcNo() )
                continue;
 
            //- the following information is sent to other processor
            //- 1. global point label in the original mesh
            //- 2. global point label in the tet mesh
            exchangeData[neiProc].append(it.key());
            exchangeData[neiProc].append(globalPointLabel[pLabel]);
        }
    }
 
    //- exchange data with other processors
    receivedData.clear();
    help::exchangeMap(exchangeData, receivedData);
 
    label counter(0);
    while( counter < receivedData.size() )
    {
        const label gpI = receivedData[counter++];
        const label tgI = receivedData[counter++];
        const label pLabel =
            nodeLabelForPoint[globalToLocalPointAddressing[gpI]];
 
        globalPointLabel[pLabel] = tgI;
    }
 
    //- set global labels for remaining points
    forAll(globalPointLabel, pI)
    {
        if( globalPointLabel[pI] == -1 )
            globalPointLabel[pI] = startPoint++;
    }
 
    //- create global to local mapping
    forAll(globalPointLabel, pI)
    {
        if( pProcs.sizeOfRow(pI) != 0 )
        {
            pAtParallelBoundaries.append(pI);
            globalToLocal.insert(globalPointLabel[pI], pI);
        }
    }
 
    //- mark vertices at parallel boundaries
    forAll(pointType_, pI)
        if( (pointType_[pI] & useType) && (pProcs.sizeOfRow(pI) != 0) )
            pointType_[pI] |= PARALLELBOUNDARY;
 
    //- create neighbour processors addressing
    if( !neiProcsPtr_ )
        neiProcsPtr_ = new DynList<label>();
    DynList<label>& neiProcs = *neiProcsPtr_;
 
    for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
        neiProcs.append(iter->first);
}
 
void partTriMesh::createBufferLayers()
{
    pointField& pts = triSurfModifier(surf_).pointsAccess();
 
    VRWGraph& pProcs = *pAtProcsPtr_;
    labelLongList& globalPointLabel = *globalPointLabelPtr_;
    Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
    const DynList<label>& neiProcs = *this->neiProcsPtr_;
 
    if( !pAtBufferLayersPtr_ )
        pAtBufferLayersPtr_ = new labelLongList();
    labelLongList& pAtBufferLayers = *pAtBufferLayersPtr_;
    pAtBufferLayers.clear();
 
    //- create the map
    std::map<label, LongList<parTriFace> > exchangeTrias;
    forAll(neiProcs, procI)
        exchangeTrias.insert
        (
            std::make_pair(neiProcs[procI], LongList<parTriFace>())
        );
 
    //- loop over triangles and add the ones having vertices at parallel
    //- boundaries for sending
    forAll(surf_, triI)
    {
        const labelledTri& pt = surf_[triI];
 
        DynList<label> sendToProcs;
        forAll(pt, i)
        {
            const label pLabel = pt[i];
 
            if( pointType_[pLabel] & PARALLELBOUNDARY )
            {
                forAllRow(pProcs, pLabel, i)
                {
                    const label neiProc = pProcs(pLabel, i);
 
                    if( neiProc == Pstream::myProcNo() )
                        continue;
 
                    sendToProcs.appendIfNotIn(neiProc);
                }
            }
        }
 
        if( sendToProcs.size() )
        {
            const parTriFace tri
            (
                globalPointLabel[pt[0]],
                globalPointLabel[pt[1]],
                globalPointLabel[pt[2]],
                triangle<point, point>(pts[pt[0]], pts[pt[1]], pts[pt[2]])
            );
 
            forAll(sendToProcs, i)
            {
                exchangeTrias[sendToProcs[i]].append(tri);
 
                forAll(pt, j)
                {
                    if( pProcs.sizeOfRow(pt[j]) == 0 )
                        pAtBufferLayers.append(pt[j]);
 
                    pProcs.appendIfNotIn(pt[j], sendToProcs[i]);
                }
            }
        }
    }
 
    //- receive triangles sent to this processor
    std::map<label, List<parTriFace> > receivedTriangles;
    help::exchangeMap(exchangeTrias, receivedTriangles);
    exchangeTrias.clear();
 
    //- add triangles into the mesh and update the addressing
    Map<label> newGlobalToLocal;
    std::map<label, point> addCoordinates;
    label nPoints = pts.size();
    for
    (
        std::map<label, List<parTriFace> >::const_iterator it =
        receivedTriangles.begin();
        it!=receivedTriangles.end();
        ++it
    )
    {
        const List<parTriFace>& receivedTrias = it->second;
 
        forAll(receivedTrias, i)
        {
            const parTriFace& tri = receivedTrias[i];
 
            DynList<label, 3> triPointLabels(3);
            for(label j=0;j<3;++j)
            {
                const label gpI = tri.globalLabelOfPoint(j);
 
                if( globalToLocal.found(gpI) )
                {
                    //- point already exists in the triangulation
                    const label pI = globalToLocal[gpI];
                    triPointLabels[j] = pI;
                }
                else if( newGlobalToLocal.found(gpI) )
                {
                    //- point is already added into the triangulation
                    triPointLabels[j] = newGlobalToLocal[gpI];
                    pProcs.appendIfNotIn(newGlobalToLocal[gpI], it->first);
                }
                else
                {
                    //- point does not exist in the triangulation
                    //- and is not yet added in
                    newGlobalToLocal.insert(gpI, nPoints);
                    triPointLabels[j] = nPoints;
 
                    point tp;
                    if( j == 0 )
                    {
                        tp = tri.trianglePoints().a();
                    }
                    else if( j == 1 )
                    {
                        tp = tri.trianglePoints().b();
                    }
                    else
                    {
                        tp = tri.trianglePoints().c();
                    }
                    addCoordinates[nPoints] = tp;
                    ++nPoints;
 
                    pointLabelInMeshSurface_.append(-1);
                    pointType_.append(NONE);
 
                    DynList<label> triAtProcs;
                    triAtProcs.append(it->first);
 
                    globalPointLabel.append(gpI);
                    triAtProcs.append(Pstream::myProcNo());
                    pProcs.appendList(triAtProcs);
                }
            }
 
            //- append tet
            surf_.appendTriangle
            (
                labelledTri
                (
                    triPointLabels[0],
                    triPointLabels[1],
                    triPointLabels[2],
                    -1
                )
            );
        }
    }
 
    //- store newly added points
    pts.setSize(nPoints);
    for
    (
        std::map<label, point>::const_iterator it=addCoordinates.begin();
        it!=addCoordinates.end();
        ++it
    )
        pts[it->first] = it->second;
 
    addCoordinates.clear();
 
    //- insert the global labels of the buffer points
    //- into the globalToLocal map
    forAllConstIter(Map<label>, newGlobalToLocal, it)
        globalToLocal.insert(it.key(), it());
 
    //- update addressing of the surface mesh
    surf_.clearAddressing();
}
 
void partTriMesh::updateBufferLayers()
{
    const pointField& points = surf_.points();
    const labelLongList& bufferLayerPoints = this->bufferLayerPoints();
    const VRWGraph& pProcs = this->pointAtProcs();
    const labelLongList& globalPointLabel = this->globalPointLabel();
    const Map<label>& globalToLocal = this->globalToLocalPointAddressing();
    const DynList<label>& neiProcs = this->neiProcs();
 
    //- create the map
    std::map<label, LongList<labelledPoint> > exchangeData;
    forAll(neiProcs, i)
        exchangeData.insert
        (
            std::make_pair(neiProcs[i], LongList<labelledPoint>())
        );
 
    //- add points into the map
    forAll(bufferLayerPoints, pI)
    {
        const label pointI = bufferLayerPoints[pI];
 
        forAllRow(pProcs, pointI, i)
        {
            const label neiProc = pProcs(pointI, i);
 
            if( neiProc == Pstream::myProcNo() )
                continue;
 
            exchangeData[neiProc].append
            (
                labelledPoint(globalPointLabel[pointI], points[pointI])
            );
        }
    }
 
    LongList<labelledPoint> receivedData;
    help::exchangeMap(exchangeData, receivedData);
 
    forAll(receivedData, i)
    {
        const labelledPoint& lp = receivedData[i];
 
        this->updateVertex
        (
            globalToLocal[lp.pointLabel()],
            lp.coordinates()
        );
    }
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// ************************************************************************* //