src/mesh/cfMesh/utilities/meshes/partTetMesh/partTetMeshParallelAddressing.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 "polyMeshGenModifier.H"
  30 #include "partTetMesh.H"
  31 #include "polyMeshGenAddressing.H"
  32 #include "helperFunctionsPar.H"
  33 #include "parPartTet.H"
  34
  35 #include <map>
  36
  37 //#define DEBUGSmooth
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 void partTetMesh::createParallelAddressing
  47 (
  48     const labelLongList& nodeLabelForPoint,
  49     const labelLongList& nodeLabelForFace,
  50     const labelLongList& nodeLabelForCell
  51 )
  52 {
  53     //- vertices marked as SMOOTH and BOUNDARY are used by the smoother
  54     const direction useType = SMOOTH + BOUNDARY;
  55
  56     //- allocate global point labels
  57     if( !globalPointLabelPtr_ )
  58         globalPointLabelPtr_ = new labelLongList();
  59     labelLongList& globalTetPointLabel = *globalPointLabelPtr_;
  60     globalTetPointLabel.setSize(points_.size());
  61     globalTetPointLabel = -1;
  62
  63     //- allocated point-processors addressing
  64     if( !pAtProcsPtr_ )
  65         pAtProcsPtr_ = new VRWGraph();
  66     VRWGraph& pProcs = *pAtProcsPtr_;
  67     pProcs.setSize(0);
  68     pProcs.setSize(points_.size());
  69
  70     //- allocate global-to-local point addressing
  71     if( !globalToLocalPointAddressingPtr_ )
  72         globalToLocalPointAddressingPtr_ = new Map<label>();
  73     Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
  74     globalToLocal.clear();
  75
  76     //- allocate storage for points at parallel boundaries
  77     if( !pAtParallelBoundariesPtr_ )
  78         pAtParallelBoundariesPtr_ = new labelLongList();
  79     labelLongList& pAtParallelBoundaries = *pAtParallelBoundariesPtr_;
  80     pAtParallelBoundaries.clear();
  81
  82     //- create point-processors addressing
  83     std::map<label, labelLongList> exchangeData;
  84     std::map<label, labelLongList>::iterator iter;
  85
  86     const polyMeshGenAddressing& addressing = origMesh_.addressingData();
  87     const Map<label>& globalToLocalPointAddressing =
  88         addressing.globalToLocalPointAddressing();
  89     const VRWGraph& pAtProcs = addressing.pointAtProcs();
  90     const DynList<label>& pNeiProcs = addressing.pointNeiProcs();
  91
  92     forAll(pNeiProcs, procI)
  93         exchangeData.insert(std::make_pair(pNeiProcs[procI], labelLongList()));
  94
  95     //- make sure that the same vertices are marked for smoothing on all procs
  96     //- this is performed by sending the labels of vertices which are not used
  97     //- for tet mesh creation and the tet mesh vertices which are not moved
  98     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
  99     {
 100         const label pI = it();
 101
 102         if(
 103             nodeLabelForPoint[pI] == -1 ||
 104             !smoothVertex_[nodeLabelForPoint[pI]]
 105         )
 106         {
 107             forAllRow(pAtProcs, pI, procI)
 108             {
 109                 const label neiProc = pAtProcs(pI, procI);
 110                 if( neiProc == Pstream::myProcNo() )
 111                     continue;
 112
 113                 exchangeData[neiProc].append(it.key());
 114             }
 115         }
 116     }
 117
 118     //- exchange data with other processors
 119     labelLongList receivedData;
 120     help::exchangeMap(exchangeData, receivedData);
 121
 122     //- set the values according to other processors
 123     forAll(receivedData, i)
 124     {
 125         const label pointI = globalToLocalPointAddressing[receivedData[i]];
 126
 127         if( nodeLabelForPoint[pointI] == -1 )
 128             continue;
 129
 130         smoothVertex_[nodeLabelForPoint[pointI]] = NONE;
 131     }
 132
 133     for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
 134         iter->second.clear();
 135
 136     //- start creating global-to-local addressing
 137     //- find the starting point labels
 138     label startPoint(0), nLocalPoints(0), nSharedPoints(0);
 139
 140     //- count the number of points at processor boundaries
 141     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
 142     {
 143         const label pI = it();
 144
 145         if( nodeLabelForPoint[pI] == -1 )
 146             continue;
 147         if( !(smoothVertex_[nodeLabelForPoint[pI]] & useType) )
 148             continue;
 149
 150         ++nSharedPoints;
 151
 152         label pMin(Pstream::myProcNo());
 153         forAllRow(pAtProcs, pI, procI)
 154             pMin = Foam::min(pMin, pAtProcs(pI, procI));
 155
 156         if( pMin == Pstream::myProcNo() )
 157             ++nLocalPoints;
 158     }
 159
 160     labelList nPointsAtProc(Pstream::nProcs());
 161     nSharedPoints -= nLocalPoints;
 162     nPointsAtProc[Pstream::myProcNo()] = points_.size() - nSharedPoints;
 163     Pstream::gatherList(nPointsAtProc);
 164     Pstream::scatterList(nPointsAtProc);
 165
 166     for(label i=0;i<Pstream::myProcNo();++i)
 167         startPoint += nPointsAtProc[i];
 168
 169     //- create global labels for points at processor boundaries
 170     forAllConstIter(Map<label>, globalToLocalPointAddressing, it)
 171     {
 172         const label pI = it();
 173
 174         if( nodeLabelForPoint[pI] == -1 )
 175             continue;
 176
 177         const label pLabel = nodeLabelForPoint[pI];
 178
 179         if( !(smoothVertex_[pLabel] & useType) )
 180             continue;
 181
 182         label pMin(Pstream::myProcNo());
 183         forAllRow(pAtProcs, pI, procI)
 184         {
 185             const label neiProc = pAtProcs(pI, procI);
 186             pProcs.append(pLabel, neiProc);
 187             pMin = Foam::min(pMin, neiProc);
 188         }
 189
 190         if( pMin != Pstream::myProcNo() )
 191             continue;
 192
 193         globalTetPointLabel[pLabel] = startPoint++;
 194
 195         forAllRow(pAtProcs, pI, procI)
 196         {
 197             const label neiProc = pAtProcs(pI, procI);
 198
 199             if( neiProc == Pstream::myProcNo() )
 200                 continue;
 201
 202             //- the following information is sent to other processor
 203             //- 1. global point label in the original mesh
 204             //- 2. global point label in the tet mesh
 205             exchangeData[neiProc].append(it.key());
 206             exchangeData[neiProc].append(globalTetPointLabel[pLabel]);
 207         }
 208     }
 209
 210     //- exchange data with other processors
 211     receivedData.clear();
 212     help::exchangeMap(exchangeData, receivedData);
 213
 214     label counter(0);
 215     while( counter < receivedData.size() )
 216     {
 217         const label gpI = receivedData[counter++];
 218         const label tgI = receivedData[counter++];
 219         const label pLabel =
 220             nodeLabelForPoint[globalToLocalPointAddressing[gpI]];
 221
 222         globalTetPointLabel[pLabel] = tgI;
 223     }
 224
 225     //- set global labels for remaining points
 226     forAll(globalTetPointLabel, pI)
 227     {
 228         if( globalTetPointLabel[pI] == -1 )
 229             globalTetPointLabel[pI] = startPoint++;
 230     }
 231
 232     //- create global to local mapping
 233     forAll(globalTetPointLabel, pI)
 234     {
 235         if( pProcs.sizeOfRow(pI) != 0 )
 236         {
 237             pAtParallelBoundaries.append(pI);
 238             globalToLocal.insert(globalTetPointLabel[pI], pI);
 239         }
 240     }
 241
 242     //- mark vertices at parallel boundaries
 243     forAll(smoothVertex_, pI)
 244         if( (smoothVertex_[pI] & useType) && (pProcs.sizeOfRow(pI) != 0) )
 245             smoothVertex_[pI] |= PARALLELBOUNDARY;
 246
 247     //- create neighbour processors addressing
 248     if( !neiProcsPtr_ )
 249         neiProcsPtr_ = new DynList<label>();
 250     DynList<label>& neiProcs = *neiProcsPtr_;
 251
 252     for(iter=exchangeData.begin();iter!=exchangeData.end();++iter)
 253         neiProcs.append(iter->first);
 254
 255     # ifdef DEBUGSmooth
 256     for(label i=0;i<Pstream::nProcs();++i)
 257     {
 258         if( i == Pstream::myProcNo() )
 259         {
 260             Pout << "globalTetPointLabel " << globalTetPointLabel << endl;
 261         }
 262
 263         returnReduce(i, sumOp<label>());
 264     }
 265
 266     returnReduce(1, sumOp<label>());
 267
 268     forAll(nodeLabelForPoint, pI)
 269     {
 270         const label tpI = nodeLabelForPoint[pI];
 271         if( tpI != -1 && globalTetPointLabel[tpI] == -1 )
 272             FatalError << "Crap1 " << tpI << abort(FatalError);
 273     }
 274
 275     returnReduce(1, sumOp<label>());
 276
 277     forAll(nodeLabelForFace, fI)
 278     {
 279         const label tpI = nodeLabelForFace[fI];
 280         if( tpI != -1 && globalTetPointLabel[tpI] == -1 )
 281         {
 282             Pout << "Face point " << tpI << " is at procs "
 283                 << pProcs[tpI] << endl;
 284             FatalError << "Crap2" << tpI << abort(FatalError);
 285         }
 286     }
 287
 288     returnReduce(1, sumOp<label>());
 289
 290     forAll(nodeLabelForCell, cI)
 291     {
 292         const label tpI = nodeLabelForCell[cI];
 293         if( tpI != -1 && globalTetPointLabel[tpI] == -1 )
 294             FatalError << "Crap3" << tpI << abort(FatalError);
 295     }
 296
 297     forAll(smoothVertex_, vI)
 298         if( smoothVertex_[vI] & partTetMesh::PARALLELBOUNDARY )
 299             Pout << "Point " << globalTetPointLabel[vI]
 300             << " is at par bnd" << endl;
 301
 302     Serr << Pstream::myProcNo() << "points " << points_ << endl;
 303     Serr << Pstream::myProcNo() << "Tets " << tets_ << endl;
 304     forAll(pProcs, pI)
 305     {
 306         if( pProcs.sizeOfRow(pI) == 0 )
 307             continue;
 308
 309         Serr << Pstream::myProcNo() << "Point " << globalTetPointLabel[pI]
 310             << " is at procs " << pProcs[pI] << " n tets "
 311             << pointTets_[pI].size() << endl;
 312     }
 313
 314     returnReduce(1, sumOp<label>());
 315     # endif
 316 }
 317
 318 void partTetMesh::createBufferLayers()
 319 {
 320     VRWGraph& pProcs = *pAtProcsPtr_;
 321     labelLongList& globalTetPointLabel = *globalPointLabelPtr_;
 322     Map<label>& globalToLocal = *globalToLocalPointAddressingPtr_;
 323     const DynList<label>& neiProcs = *this->neiProcsPtr_;
 324
 325     if( !pAtBufferLayersPtr_ )
 326         pAtBufferLayersPtr_ = new labelLongList();
 327     labelLongList& pAtBufferLayers = *pAtBufferLayersPtr_;
 328     pAtBufferLayers.clear();
 329
 330     //- create the map
 331     std::map<label, LongList<parPartTet> > exchangeTets;
 332     forAll(neiProcs, procI)
 333         exchangeTets.insert
 334         (
 335             std::make_pair(neiProcs[procI], LongList<parPartTet>())
 336         );
 337
 338     //- go through the tets and add the ones having vertices at parallel
 339     //- boundaries for sending
 340     forAll(tets_, tetI)
 341     {
 342         const partTet& pt = tets_[tetI];
 343
 344         DynList<label> sendToProcs;
 345         forAll(pt, i)
 346         {
 347             const label pLabel = pt[i];
 348
 349             if( smoothVertex_[pLabel] & PARALLELBOUNDARY )
 350             {
 351                 forAllRow(pProcs, pLabel, i)
 352                 {
 353                     const label neiProc = pProcs(pLabel, i);
 354
 355                     if( neiProc == Pstream::myProcNo() )
 356                         continue;
 357
 358                     sendToProcs.appendIfNotIn(neiProc);
 359                 }
 360             }
 361         }
 362
 363         if( sendToProcs.size() )
 364         {
 365             const parPartTet tet
 366             (
 367                 labelledPoint(globalTetPointLabel[pt[0]], points_[pt[0]]),
 368                 labelledPoint(globalTetPointLabel[pt[1]], points_[pt[1]]),
 369                 labelledPoint(globalTetPointLabel[pt[2]], points_[pt[2]]),
 370                 labelledPoint(globalTetPointLabel[pt[3]], points_[pt[3]])
 371             );
 372
 373             forAll(sendToProcs, i)
 374             {
 375                 exchangeTets[sendToProcs[i]].append(tet);
 376
 377                 forAll(pt, j)
 378                 {
 379                     if( pProcs.sizeOfRow(pt[j]) == 0 )
 380                         pAtBufferLayers.append(pt[j]);
 381
 382                     pProcs.appendIfNotIn(pt[j], sendToProcs[i]);
 383                 }
 384             }
 385         }
 386     }
 387
 388     LongList<parPartTet> receivedTets;
 389     help::exchangeMap(exchangeTets, receivedTets);
 390     exchangeTets.clear();
 391
 392     Map<label> newGlobalToLocal;
 393     forAll(receivedTets, i)
 394     {
 395         const parPartTet& tet = receivedTets[i];
 396
 397         DynList<label> tetPointLabels;
 398         for(label j=0;j<4;++j)
 399         {
 400             const label gpI = tet[j].pointLabel();
 401
 402             if( globalToLocal.found(gpI) )
 403             {
 404                 const label pI = globalToLocal[gpI];
 405                 pointTets_.append(pI, tets_.size());
 406                 tetPointLabels.append(pI);
 407             }
 408             else if( newGlobalToLocal.found(gpI) )
 409             {
 410                 tetPointLabels.append(newGlobalToLocal[gpI]);
 411             }
 412             else
 413             {
 414                 newGlobalToLocal.insert(gpI, points_.size());
 415                 tetPointLabels.append(points_.size());
 416                 points_.append(tet[j].coordinates());
 417                 nodeLabelInOrigMesh_.append(-1);
 418                 smoothVertex_.append(NONE);
 419                 DynList<label> helper;
 420                 helper.append(tets_.size());
 421                 pointTets_.appendList(helper);
 422
 423                 globalTetPointLabel.append(gpI);
 424                 helper[0] = Pstream::myProcNo();
 425                 pProcs.appendList(helper);
 426             }
 427         }
 428
 429         //- append tet
 430         tets_.append
 431         (
 432             partTet
 433             (
 434                 tetPointLabels[0],
 435                 tetPointLabels[1],
 436                 tetPointLabels[2],
 437                 tetPointLabels[3]
 438             )
 439         );
 440     }
 441
 442     //- insert the global labels of the buffer points
 443     //- into the globalToLocal map
 444     forAllConstIter(Map<label>, newGlobalToLocal, it)
 445         globalToLocal.insert(it.key(), it());
 446
 447     # ifdef DEBUGSmooth
 448     for(label i=0;i<Pstream::nProcs();++i)
 449     {
 450         if( Pstream::myProcNo() == i )
 451         {
 452             forAllConstIter(Map<label>, globalToLocal, it)
 453             {
 454                 DynList<label> np;
 455                 if( it() < pProcs.size() )
 456                 forAllRow(pProcs, it(), j)
 457                     np.append(pProcs(it(), j));
 458
 459                 Pout << "Tet mesh point " << it() << " has global label "
 460                     << it.key() << " and is located at procs "
 461                     << np << endl;
 462             }
 463         }
 464
 465         returnReduce(1, sumOp<label>());
 466     }
 467     # endif
 468 }
 469
 470 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 471
 472 } // End namespace Foam
 473
 474 // ************************************************************************* //