src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshParallelAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "demandDrivenData.H"
  29 #include "meshSurfacePartitioner.H"
  30 #include "partTriMesh.H"
  31 #include "triSurfModifier.H"
  32 #include "meshSurfaceEngine.H"
  33 #include "helperFunctionsPar.H"
  34 #include "parTriFace.H"
  35
  36 #include <map>
  37
  38 //#define DEBUGSmooth
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 void partTriMesh::createParallelAddressing
  48 (
  49     const labelList& nodeLabelForPoint,
  50     const labelList& nodeLabelForFace
  51 )
  52 {
  53     const meshSurfaceEngine& mse = mPart_.surfaceEngine();
  54
  55     const pointField& pts = surf_.points();
  56
  57     //- vertices marked as SMOOTH are used by the smoother
  58     const direction useType = SMOOTH;
  59
  60     //- allocate global point labels
  61     if( !globalPointLabelPtr_ )
  62         globalPointLabelPtr_ = new labelLongList();
  63     labelLongList& globalPointLabel = *globalPointLabelPtr_;
  64     globalPointLabel.setSize(pts.size());
  65     globalPointLabel = -1;
  66
  67     //- allocated point-processors addressing
  68     if( !pAtProcsPtr_ )
  69         pAtProcsPtr_ = new VRWGraph();
  70     VRWGraph& pProcs = *pAtProcsPtr_;
  71     pProcs.setSize(0);
  72     pProcs.setSize(pts.size());
  73
  74     //- allocate global-to-local point addressing
  75     if( !globalToLocalPointAddressingPtr_ )
  76         globalToLocalPointAddressingPtr_ = new Map<label>();
  77     Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
  78     globalToLocal.clear();
  79
  80     //- allocate storage for points at parallel boundaries
  81     if( !pAtParallelBoundariesPtr_ )
  82         pAtParallelBoundariesPtr_ = new labelLongList();
  83     labelLongList& pAtParallelBoundaries = *pAtParallelBoundariesPtr_;
  84     pAtParallelBoundaries.clear();
  85
  86     //- create point-processors addressing
  87     std::map<label, labelLongList> exchangeData;
  88     std::map<label, labelLongList>::iterator iter;
  89
  90     const Map<label>& globalToLocalPointAddressing =
  91         mse.globalToLocalBndPointAddressing();
  92     const VRWGraph& pAtProcs = mse.bpAtProcs();
  93     const DynList<label>& pNeiProcs = mse.bpNeiProcs();
  94
  95     forAll(pNeiProcs, procI)
  96         exchangeData.insert(std::make_pair(pNeiProcs[procI], labelLongList()));
  97
  98     //- make sure that the same vertices are marked for smoothing on all procs
  99     //- this is performed by sending the labels of vertices which are not used
 100     //- for tet mesh creation and the tet mesh vertices which are not moved
 101     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
 102     {
 103         const label pI = it();
 104
 105         if(
 106             nodeLabelForPoint[pI] == -1 ||
 107             !pointType_[nodeLabelForPoint[pI]]
 108         )
 109         {
 110             forAllRow(pAtProcs, pI, procI)
 111             {
 112                 const label neiProc = pAtProcs(pI, procI);
 113                 if( neiProc == Pstream::myProcNo() )
 114                     continue;
 115
 116                 exchangeData[neiProc].append(it.key());
 117             }
 118         }
 119     }
 120
 121     //- exchange data with other processors
 122     labelLongList receivedData;
 123     help::exchangeMap(exchangeData, receivedData);
 124
 125     //- set the values according to other processors
 126     forAll(receivedData, i)
 127     {
 128         const label pointI = globalToLocalPointAddressing[receivedData[i]];
 129
 130         if( nodeLabelForPoint[pointI] == -1 )
 131             continue;
 132
 133         pointType_[nodeLabelForPoint[pointI]] = NONE;
 134     }
 135
 136     for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
 137         iter->second.clear();
 138
 139     //- start creating global-to-local addressing
 140     //- find the starting point labels
 141     label startPoint(0), nLocalPoints(0), nSharedPoints(0);
 142
 143     //- count the number of points at processor boundaries
 144     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
 145     {
 146         const label pI = it();
 147
 148         if( nodeLabelForPoint[pI] == -1 )
 149             continue;
 150         if( !(pointType_[nodeLabelForPoint[pI]] & useType) )
 151             continue;
 152
 153         ++nSharedPoints;
 154
 155         label pMin(Pstream::myProcNo());
 156         forAllRow(pAtProcs, pI, procI)
 157             pMin = Foam::min(pMin, pAtProcs(pI, procI));
 158
 159         if( pMin == Pstream::myProcNo() )
 160             ++nLocalPoints;
 161     }
 162
 163     labelList nPointsAtProc(Pstream::nProcs());
 164     nSharedPoints -= nLocalPoints;
 165     nPointsAtProc[Pstream::myProcNo()] = pts.size() - nSharedPoints;
 166     Pstream::gatherList(nPointsAtProc);
 167     Pstream::scatterList(nPointsAtProc);
 168
 169     for(label i=0;i<Pstream::myProcNo();++i)
 170         startPoint += nPointsAtProc[i];
 171
 172     //- create global labels for points at processor boundaries
 173     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
 174     {
 175         const label pI = it();
 176
 177         if( nodeLabelForPoint[pI] == -1 )
 178             continue;
 179
 180         const label pLabel = nodeLabelForPoint[pI];
 181
 182         if( !(pointType_[pLabel] & useType) )
 183             continue;
 184
 185         label pMin(Pstream::myProcNo());
 186         forAllRow(pAtProcs, pI, procI)
 187         {
 188             const label neiProc = pAtProcs(pI, procI);
 189             pProcs.append(pLabel, neiProc);
 190             pMin = Foam::min(pMin, neiProc);
 191         }
 192
 193         if( pMin != Pstream::myProcNo() )
 194             continue;
 195
 196         globalPointLabel[pLabel] = startPoint++;
 197
 198         forAllRow(pAtProcs, pI, procI)
 199         {
 200             const label neiProc = pAtProcs(pI, procI);
 201
 202             if( neiProc == Pstream::myProcNo() )
 203                 continue;
 204
 205             //- the following information is sent to other processor
 206             //- 1. global point label in the original mesh
 207             //- 2. global point label in the tet mesh
 208             exchangeData[neiProc].append(it.key());
 209             exchangeData[neiProc].append(globalPointLabel[pLabel]);
 210         }
 211     }
 212
 213     //- exchange data with other processors
 214     receivedData.clear();
 215     help::exchangeMap(exchangeData, receivedData);
 216
 217     label counter(0);
 218     while( counter < receivedData.size() )
 219     {
 220         const label gpI = receivedData[counter++];
 221         const label tgI = receivedData[counter++];
 222         const label pLabel =
 223             nodeLabelForPoint[globalToLocalPointAddressing[gpI]];
 224
 225         globalPointLabel[pLabel] = tgI;
 226     }
 227
 228     //- set global labels for remaining points
 229     forAll(globalPointLabel, pI)
 230     {
 231         if( globalPointLabel[pI] == -1 )
 232             globalPointLabel[pI] = startPoint++;
 233     }
 234
 235     //- create global to local mapping
 236     forAll(globalPointLabel, pI)
 237     {
 238         if( pProcs.sizeOfRow(pI) != 0 )
 239         {
 240             pAtParallelBoundaries.append(pI);
 241             globalToLocal.insert(globalPointLabel[pI], pI);
 242         }
 243     }
 244
 245     //- mark vertices at parallel boundaries
 246     forAll(pointType_, pI)
 247         if( (pointType_[pI] & useType) && (pProcs.sizeOfRow(pI) != 0) )
 248             pointType_[pI] |= PARALLELBOUNDARY;
 249
 250     //- create neighbour processors addressing
 251     if( !neiProcsPtr_ )
 252         neiProcsPtr_ = new DynList<label>();
 253     DynList<label>& neiProcs = *neiProcsPtr_;
 254
 255     for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
 256         neiProcs.append(iter->first);
 257 }
 258
 259 void partTriMesh::createBufferLayers()
 260 {
 261     pointField& pts = triSurfModifier(surf_).pointsAccess();
 262
 263     VRWGraph& pProcs = *pAtProcsPtr_;
 264     labelLongList& globalPointLabel = *globalPointLabelPtr_;
 265     Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
 266     const DynList<label>& neiProcs = *this->neiProcsPtr_;
 267
 268     if( !pAtBufferLayersPtr_ )
 269         pAtBufferLayersPtr_ = new labelLongList();
 270     labelLongList& pAtBufferLayers = *pAtBufferLayersPtr_;
 271     pAtBufferLayers.clear();
 272
 273     //- create the map
 274     std::map<label, LongList<parTriFace> > exchangeTrias;
 275     forAll(neiProcs, procI)
 276         exchangeTrias.insert
 277         (
 278             std::make_pair(neiProcs[procI], LongList<parTriFace>())
 279         );
 280
 281     //- loop over triangles and add the ones having vertices at parallel
 282     //- boundaries for sending
 283     forAll(surf_, triI)
 284     {
 285         const labelledTri& pt = surf_[triI];
 286
 287         DynList<label> sendToProcs;
 288         forAll(pt, i)
 289         {
 290             const label pLabel = pt[i];
 291
 292             if( pointType_[pLabel] & PARALLELBOUNDARY )
 293             {
 294                 forAllRow(pProcs, pLabel, i)
 295                 {
 296                     const label neiProc = pProcs(pLabel, i);
 297
 298                     if( neiProc == Pstream::myProcNo() )
 299                         continue;
 300
 301                     sendToProcs.appendIfNotIn(neiProc);
 302                 }
 303             }
 304         }
 305
 306         if( sendToProcs.size() )
 307         {
 308             const parTriFace tri
 309             (
 310                 globalPointLabel[pt[0]],
 311                 globalPointLabel[pt[1]],
 312                 globalPointLabel[pt[2]],
 313                 triangle<point, point>(pts[pt[0]], pts[pt[1]], pts[pt[2]])
 314             );
 315
 316             forAll(sendToProcs, i)
 317             {
 318                 exchangeTrias[sendToProcs[i]].append(tri);
 319
 320                 forAll(pt, j)
 321                 {
 322                     if( pProcs.sizeOfRow(pt[j]) == 0 )
 323                         pAtBufferLayers.append(pt[j]);
 324
 325                     pProcs.appendIfNotIn(pt[j], sendToProcs[i]);
 326                 }
 327             }
 328         }
 329     }
 330
 331     //- receive triangles sent to this processor
 332     std::map<label, List<parTriFace> > receivedTriangles;
 333     help::exchangeMap(exchangeTrias, receivedTriangles);
 334     exchangeTrias.clear();
 335
 336     //- add triangles into the mesh and update the addressing
 337     Map<label> newGlobalToLocal;
 338     std::map<label, point> addCoordinates;
 339     label nPoints = pts.size();
 340     for
 341     (
 342         std::map<label, List<parTriFace> >::const_iterator it =
 343         receivedTriangles.begin();
 344         it!=receivedTriangles.end();
 345         ++it
 346     )
 347     {
 348         const List<parTriFace>& receivedTrias = it->second;
 349
 350         forAll(receivedTrias, i)
 351         {
 352             const parTriFace& tri = receivedTrias[i];
 353
 354             DynList<label, 3> triPointLabels(3);
 355             for(label j=0;j<3;++j)
 356             {
 357                 const label gpI = tri.globalLabelOfPoint(j);
 358
 359                 if( globalToLocal.found(gpI) )
 360                 {
 361                     //- point already exists in the triangulation
 362                     const label pI = globalToLocal[gpI];
 363                     triPointLabels[j] = pI;
 364                 }
 365                 else if( newGlobalToLocal.found(gpI) )
 366                 {
 367                     //- point is already added into the triangulation
 368                     triPointLabels[j] = newGlobalToLocal[gpI];
 369                     pProcs.appendIfNotIn(newGlobalToLocal[gpI], it->first);
 370                 }
 371                 else
 372                 {
 373                     //- point does not exist in the triangulation
 374                     //- and is not yet added in
 375                     newGlobalToLocal.insert(gpI, nPoints);
 376                     triPointLabels[j] = nPoints;
 377
 378                     point tp;
 379                     if( j == 0 )
 380                     {
 381                         tp = tri.trianglePoints().a();
 382                     }
 383                     else if( j == 1 )
 384                     {
 385                         tp = tri.trianglePoints().b();
 386                     }
 387                     else
 388                     {
 389                         tp = tri.trianglePoints().c();
 390                     }
 391                     addCoordinates[nPoints] = tp;
 392                     ++nPoints;
 393
 394                     pointLabelInMeshSurface_.append(-1);
 395                     pointType_.append(NONE);
 396
 397                     DynList<label> triAtProcs;
 398                     triAtProcs.append(it->first);
 399
 400                     globalPointLabel.append(gpI);
 401                     triAtProcs.append(Pstream::myProcNo());
 402                     pProcs.appendList(triAtProcs);
 403                 }
 404             }
 405
 406             //- append tet
 407             surf_.appendTriangle
 408             (
 409                 labelledTri
 410                 (
 411                     triPointLabels[0],
 412                     triPointLabels[1],
 413                     triPointLabels[2],
 414                     -1
 415                 )
 416             );
 417         }
 418     }
 419
 420     //- store newly added points
 421     pts.setSize(nPoints);
 422     for
 423     (
 424         std::map<label, point>::const_iterator it=addCoordinates.begin();
 425         it!=addCoordinates.end();
 426         ++it
 427     )
 428         pts[it->first] = it->second;
 429
 430     addCoordinates.clear();
 431
 432     //- insert the global labels of the buffer points
 433     //- into the globalToLocal map
 434     forAllConstIter(Map<label>, newGlobalToLocal, it)
 435         globalToLocal.insert(it.key(), it());
 436
 437     //- update addressing of the surface mesh
 438     surf_.clearAddressing();
 439 }
 440
 441 void partTriMesh::updateBufferLayers()
 442 {
 443     const pointField& points = surf_.points();
 444     const labelLongList& bufferLayerPoints = this->bufferLayerPoints();
 445     const VRWGraph& pProcs = this->pointAtProcs();
 446     const labelLongList& globalPointLabel = this->globalPointLabel();
 447     const Map<label>& globalToLocal = this->globalToLocalPointAddressing();
 448     const DynList<label>& neiProcs = this->neiProcs();
 449
 450     //- create the map
 451     std::map<label, LongList<labelledPoint> > exchangeData;
 452     forAll(neiProcs, i)
 453         exchangeData.insert
 454         (
 455             std::make_pair(neiProcs[i], LongList<labelledPoint>())
 456         );
 457
 458     //- add points into the map
 459     forAll(bufferLayerPoints, pI)
 460     {
 461         const label pointI = bufferLayerPoints[pI];
 462
 463         forAllRow(pProcs, pointI, i)
 464         {
 465             const label neiProc = pProcs(pointI, i);
 466
 467             if( neiProc == Pstream::myProcNo() )
 468                 continue;
 469
 470             exchangeData[neiProc].append
 471             (
 472                 labelledPoint(globalPointLabel[pointI], points[pointI])
 473             );
 474         }
 475     }
 476
 477     LongList<labelledPoint> receivedData;
 478     help::exchangeMap(exchangeData, receivedData);
 479
 480     forAll(receivedData, i)
 481     {
 482         const labelledPoint& lp = receivedData[i];
 483
 484         this->updateVertex
 485         (
 486             globalToLocal[lp.pointLabel()],
 487             lp.coordinates()
 488         );
 489     }
 490 }
 491
 492 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 493
 494 } // End namespace Foam
 495
 496 // ************************************************************************* //