applications/utilities/mesh/generation/blockMesh/createMergeList.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM 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 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "blockMesh.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 Foam::labelList Foam::blockMesh::createMergeList()
  32 {
  33     Info<< nl << "Creating merge list " << flush;
  34
  35     labelList MergeList(nPoints_, -1);
  36
  37     blockMesh& blocks = *this;
  38
  39     const pointField& blockPoints = topology().points();
  40     const cellList& blockCells = topology().cells();
  41     const faceList& blockFaces = topology().faces();
  42     const labelList& faceOwnerBlocks = topology().faceOwner();
  43
  44     // For efficiency, create merge pairs in the first pass
  45     labelListListList glueMergePairs(blockFaces.size());
  46
  47     const labelList& faceNeighbourBlocks = topology().faceNeighbour();
  48
  49     forAll(blockFaces, blockFaceLabel)
  50     {
  51         label blockPlabel = faceOwnerBlocks[blockFaceLabel];
  52         const pointField& blockPpoints = blocks[blockPlabel].points();
  53         const labelList& blockPfaces = blockCells[blockPlabel];
  54
  55         bool foundFace = false;
  56         label blockPfaceLabel;
  57         for
  58         (
  59             blockPfaceLabel = 0;
  60             blockPfaceLabel < blockPfaces.size();
  61             blockPfaceLabel++
  62         )
  63         {
  64             if
  65             (
  66                 blockFaces[blockPfaces[blockPfaceLabel]]
  67              == blockFaces[blockFaceLabel]
  68             )
  69             {
  70                 foundFace = true;
  71                 break;
  72             }
  73         }
  74
  75         if (!foundFace)
  76         {
  77             FatalErrorIn("blockMesh::createMergeList()")
  78                 << "Cannot find merge face for block " << blockPlabel
  79                 << exit(FatalError);
  80         };
  81
  82         const labelListList& blockPfaceFaces =
  83             blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
  84
  85         labelListList& curPairs = glueMergePairs[blockFaceLabel];
  86         curPairs.setSize(blockPfaceFaces.size());
  87
  88         // Calculate sqr of the merge tolerance as 1/10th of the min sqr
  89         // point to point distance on the block face.
  90         // At the same time merge collated points on the block's faces
  91         // (removes boundary poles etc.)
  92         // Collated points detected by initally taking a constant factor of
  93         // the size of the block.
  94
  95         boundBox bb(blockCells[blockPlabel].points(blockFaces, blockPoints));
  96         const scalar mergeSqrDist = SMALL*magSqr(bb.span());
  97
  98         // This is an N^2 algorithm
  99
 100         scalar sqrMergeTol = GREAT;
 101
 102         forAll(blockPfaceFaces, blockPfaceFaceLabel)
 103         {
 104             const labelList& blockPfaceFacePoints
 105                 = blockPfaceFaces[blockPfaceFaceLabel];
 106
 107             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
 108             {
 109                 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel2)
 110                 {
 111                     if (blockPfaceFacePointLabel != blockPfaceFacePointLabel2)
 112                     {
 113                         scalar magSqrDist = magSqr
 114                         (
 115                             blockPpoints[blockPfaceFacePoints
 116                             [blockPfaceFacePointLabel]]
 117                             - blockPpoints[blockPfaceFacePoints
 118                             [blockPfaceFacePointLabel2]]
 119                         );
 120
 121                         if (magSqrDist < mergeSqrDist)
 122                         {
 123                             label PpointLabel =
 124                                 blockPfaceFacePoints[blockPfaceFacePointLabel]
 125                               + blockOffsets_[blockPlabel];
 126
 127                             label PpointLabel2 =
 128                                 blockPfaceFacePoints[blockPfaceFacePointLabel2]
 129                               + blockOffsets_[blockPlabel];
 130
 131                             label minPP2 = min(PpointLabel, PpointLabel2);
 132
 133                             if (MergeList[PpointLabel] != -1)
 134                             {
 135                                 minPP2 = min(minPP2, MergeList[PpointLabel]);
 136                             }
 137
 138                             if (MergeList[PpointLabel2] != -1)
 139                             {
 140                                 minPP2 = min(minPP2, MergeList[PpointLabel2]);
 141                             }
 142
 143                             MergeList[PpointLabel] = MergeList[PpointLabel2]
 144                                 = minPP2;
 145                         }
 146                         else
 147                         {
 148                             sqrMergeTol = min(sqrMergeTol, magSqrDist);
 149                         }
 150                     }
 151                 }
 152             }
 153         }
 154
 155         sqrMergeTol /= 10.0;
 156
 157
 158         if (topology().isInternalFace(blockFaceLabel))
 159         {
 160         label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
 161         const pointField& blockNpoints = blocks[blockNlabel].points();
 162         const labelList& blockNfaces = blockCells[blockNlabel];
 163
 164         foundFace = false;
 165         label blockNfaceLabel;
 166         for
 167         (
 168             blockNfaceLabel = 0;
 169             blockNfaceLabel < blockNfaces.size();
 170             blockNfaceLabel++
 171         )
 172         {
 173             if
 174             (
 175                 blockFaces[blockNfaces[blockNfaceLabel]]
 176              == blockFaces[blockFaceLabel]
 177             )
 178             {
 179                 foundFace = true;
 180                 break;
 181             }
 182         }
 183
 184         if (!foundFace)
 185         {
 186             FatalErrorIn("blockMesh::createMergeList()")
 187                 << "Cannot find merge face for block " << blockNlabel
 188                 << exit(FatalError);
 189         };
 190
 191         const labelListList& blockNfaceFaces =
 192             blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
 193
 194         if (blockPfaceFaces.size() != blockNfaceFaces.size())
 195         {
 196             FatalErrorIn("blockMesh::createMergeList()")
 197                 << "Inconsistent number of faces between block pair "
 198                 << blockPlabel << " and " << blockNlabel
 199                 << exit(FatalError);
 200         }
 201
 202
 203         bool found = false;
 204
 205         // N-squared point search over all points of all faces of
 206         // master block over all point of all faces of slave block
 207         forAll(blockPfaceFaces, blockPfaceFaceLabel)
 208         {
 209             const labelList& blockPfaceFacePoints
 210                 = blockPfaceFaces[blockPfaceFaceLabel];
 211
 212             labelList& cp = curPairs[blockPfaceFaceLabel];
 213             cp.setSize(blockPfaceFacePoints.size());
 214             cp = -1;
 215
 216             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
 217             {
 218                 found = false;
 219
 220                 forAll(blockNfaceFaces, blockNfaceFaceLabel)
 221                 {
 222                     const labelList& blockNfaceFacePoints
 223                         = blockNfaceFaces[blockNfaceFaceLabel];
 224
 225                     forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
 226                     {
 227                         if
 228                         (
 229                             magSqr
 230                             (
 231                                 blockPpoints
 232                                 [blockPfaceFacePoints[blockPfaceFacePointLabel]]
 233                               - blockNpoints
 234                                 [blockNfaceFacePoints[blockNfaceFacePointLabel]]
 235                             ) < sqrMergeTol
 236                         )
 237                         {
 238                             // Found a new pair
 239                             found = true;
 240
 241                             cp[blockPfaceFacePointLabel] =
 242                                 blockNfaceFacePoints[blockNfaceFacePointLabel];
 243
 244                             label PpointLabel =
 245                                 blockPfaceFacePoints[blockPfaceFacePointLabel]
 246                               + blockOffsets_[blockPlabel];
 247
 248                             label NpointLabel =
 249                                 blockNfaceFacePoints[blockNfaceFacePointLabel]
 250                               + blockOffsets_[blockNlabel];
 251
 252                             label minPN = min(PpointLabel, NpointLabel);
 253
 254                             if (MergeList[PpointLabel] != -1)
 255                             {
 256                                 minPN = min(minPN, MergeList[PpointLabel]);
 257                             }
 258
 259                             if (MergeList[NpointLabel] != -1)
 260                             {
 261                                 minPN = min(minPN, MergeList[NpointLabel]);
 262                             }
 263
 264                             MergeList[PpointLabel] = MergeList[NpointLabel]
 265                                 = minPN;
 266                         }
 267                     }
 268                 }
 269             }
 270             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
 271             {
 272                 if (cp[blockPfaceFacePointLabel] == -1)
 273                 {
 274                     FatalErrorIn("blockMesh::createMergeList()")
 275                         << "Inconsistent point locations between block pair "
 276                         << blockPlabel << " and " << blockNlabel << nl
 277                         << "    probably due to inconsistent grading."
 278                         << exit(FatalError);
 279                 }
 280             }
 281         }
 282         }
 283     }
 284
 285
 286     const faceList::subList blockInternalFaces
 287     (
 288         blockFaces,
 289         topology().nInternalFaces()
 290     );
 291
 292     bool changedPointMerge = false;
 293     label nPasses = 0;
 294
 295     do
 296     {
 297         changedPointMerge = false;
 298         nPasses++;
 299
 300         forAll(blockInternalFaces, blockFaceLabel)
 301         {
 302             label blockPlabel = faceOwnerBlocks[blockFaceLabel];
 303             label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
 304
 305             const labelList& blockPfaces = blockCells[blockPlabel];
 306             const labelList& blockNfaces = blockCells[blockNlabel];
 307
 308             const labelListList& curPairs = glueMergePairs[blockFaceLabel];
 309
 310             bool foundFace = false;
 311             label blockPfaceLabel;
 312             for
 313             (
 314                 blockPfaceLabel = 0;
 315                 blockPfaceLabel < blockPfaces.size();
 316                 blockPfaceLabel++
 317             )
 318             {
 319                 if
 320                 (
 321                     blockFaces[blockPfaces[blockPfaceLabel]]
 322                  == blockInternalFaces[blockFaceLabel]
 323                 )
 324                 {
 325                     foundFace = true;
 326                     break;
 327                 }
 328             }
 329
 330             foundFace = false;
 331             label blockNfaceLabel;
 332             for
 333             (
 334                 blockNfaceLabel = 0;
 335                 blockNfaceLabel < blockNfaces.size();
 336                 blockNfaceLabel++
 337             )
 338             {
 339                 if
 340                 (
 341                     blockFaces[blockNfaces[blockNfaceLabel]]
 342                  == blockInternalFaces[blockFaceLabel]
 343                 )
 344                 {
 345                     foundFace = true;
 346                     break;
 347                 }
 348             }
 349
 350             const labelListList& blockPfaceFaces =
 351                 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
 352
 353             forAll(blockPfaceFaces, blockPfaceFaceLabel)
 354             {
 355                 const labelList& blockPfaceFacePoints
 356                     = blockPfaceFaces[blockPfaceFaceLabel];
 357
 358                 const labelList& cp = curPairs[blockPfaceFaceLabel];
 359
 360                 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
 361                 {
 362                     label PpointLabel =
 363                         blockPfaceFacePoints[blockPfaceFacePointLabel]
 364                       + blockOffsets_[blockPlabel];
 365
 366                     label NpointLabel =
 367                         cp[blockPfaceFacePointLabel]
 368                       + blockOffsets_[blockNlabel];
 369
 370                     if
 371                     (
 372                         MergeList[PpointLabel]
 373                      != MergeList[NpointLabel]
 374                     )
 375                     {
 376                         changedPointMerge = true;
 377
 378                         MergeList[PpointLabel]
 379                       = MergeList[NpointLabel]
 380                       = min
 381                         (
 382                             MergeList[PpointLabel],
 383                             MergeList[NpointLabel]
 384                         );
 385                     }
 386                 }
 387             }
 388         }
 389         Info << "." << flush;
 390
 391         if (nPasses > 100)
 392         {
 393             FatalErrorIn("blockMesh::createMergeList()")
 394                 << "Point merging failed after max number of passes."
 395                 << abort(FatalError);
 396         }
 397     }
 398     while (changedPointMerge);
 399     Info << endl;
 400
 401     forAll(blockInternalFaces, blockFaceLabel)
 402     {
 403         label blockPlabel = faceOwnerBlocks[blockFaceLabel];
 404         label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
 405
 406         const labelList& blockPfaces = blockCells[blockPlabel];
 407         const labelList& blockNfaces = blockCells[blockNlabel];
 408
 409         const pointField& blockPpoints = blocks[blockPlabel].points();
 410         const pointField& blockNpoints = blocks[blockNlabel].points();
 411
 412         bool foundFace = false;
 413         label blockPfaceLabel;
 414         for
 415         (
 416             blockPfaceLabel = 0;
 417             blockPfaceLabel < blockPfaces.size();
 418             blockPfaceLabel++
 419         )
 420         {
 421             if
 422             (
 423                 blockFaces[blockPfaces[blockPfaceLabel]]
 424              == blockInternalFaces[blockFaceLabel]
 425             )
 426             {
 427                 foundFace = true;
 428                 break;
 429             }
 430         }
 431
 432         if (!foundFace)
 433         {
 434             FatalErrorIn("blockMesh::createMergeList()")
 435                 << "Cannot find merge face for block " << blockPlabel
 436                 << exit(FatalError);
 437         };
 438
 439         foundFace = false;
 440         label blockNfaceLabel;
 441         for
 442         (
 443             blockNfaceLabel = 0;
 444             blockNfaceLabel < blockNfaces.size();
 445             blockNfaceLabel++
 446         )
 447         {
 448             if
 449             (
 450                 blockFaces[blockNfaces[blockNfaceLabel]]
 451              == blockInternalFaces[blockFaceLabel]
 452             )
 453             {
 454                 foundFace = true;
 455                 break;
 456             }
 457         }
 458
 459         if (!foundFace)
 460         {
 461             FatalErrorIn("blockMesh::createMergeList()")
 462                 << "Cannot find merge face for block " << blockNlabel
 463                 << exit(FatalError);
 464         };
 465
 466         const labelListList& blockPfaceFaces =
 467             blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
 468
 469         const labelListList& blockNfaceFaces =
 470             blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
 471
 472         forAll(blockPfaceFaces, blockPfaceFaceLabel)
 473         {
 474             const labelList& blockPfaceFacePoints
 475                 = blockPfaceFaces[blockPfaceFaceLabel];
 476
 477             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
 478             {
 479                 label PpointLabel =
 480                     blockPfaceFacePoints[blockPfaceFacePointLabel]
 481                   + blockOffsets_[blockPlabel];
 482
 483                 if (MergeList[PpointLabel] == -1)
 484                 {
 485                     FatalErrorIn("blockMesh::createMergeList()")
 486                         << "Unable to merge point "
 487                         << blockPfaceFacePointLabel
 488                         << ' ' << blockPpoints[blockPfaceFacePointLabel]
 489                         << " of face "
 490                         << blockPfaceLabel
 491                         << " of block "
 492                         << blockPlabel
 493                         << exit(FatalError);
 494                 }
 495             }
 496         }
 497
 498         forAll(blockNfaceFaces, blockNfaceFaceLabel)
 499         {
 500             const labelList& blockNfaceFacePoints
 501                 = blockNfaceFaces[blockNfaceFaceLabel];
 502
 503             forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
 504             {
 505                 label NpointLabel =
 506                     blockNfaceFacePoints[blockNfaceFacePointLabel]
 507                   + blockOffsets_[blockNlabel];
 508
 509                 if (MergeList[NpointLabel] == -1)
 510                 {
 511                     FatalErrorIn("blockMesh::createMergeList()")
 512                         << "unable to merge point "
 513                         << blockNfaceFacePointLabel
 514                         << ' ' << blockNpoints[blockNfaceFacePointLabel]
 515                         << " of face "
 516                         << blockNfaceLabel
 517                         << " of block "
 518                         << blockNlabel
 519                         << exit(FatalError);
 520                 }
 521             }
 522         }
 523     }
 524
 525
 526     // sort merge list to return new point label (in new shorter list)
 527     // given old point label
 528     label newPointLabel = 0;
 529
 530     forAll(MergeList, pointLabel)
 531     {
 532         if (MergeList[pointLabel] > pointLabel)
 533         {
 534             FatalErrorIn("blockMesh::createMergeList()")
 535                 << "ouch" << exit(FatalError);
 536         }
 537
 538         if
 539         (
 540             (MergeList[pointLabel] == -1)
 541          || MergeList[pointLabel] == pointLabel
 542         )
 543         {
 544             MergeList[pointLabel] = newPointLabel;
 545             newPointLabel++;
 546         }
 547         else
 548         {
 549             MergeList[pointLabel] = MergeList[MergeList[pointLabel]];
 550         }
 551     }
 552
 553     nPoints_ = newPointLabel;
 554
 555
 556     return MergeList;
 557 }
 558
 559 // ************************************************************************* //