src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Calculating cut faces of the enriched patch, together with the addressing
  26     into master and slave patch.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "enrichedPatch.H"
  31 #include "boolList.H"
  32 #include "DynamicList.H"
  33 #include "labelPair.H"
  34 #include "primitiveMesh.H"
  35 #include "HashSet.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 // If the cut face gets more than this number of points, it will be checked
  40 const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
  41
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 // Index of debug signs:
  46 // x - skip a point
  47 // l - left turn
  48 // r - right turn
  49
  50 void Foam::enrichedPatch::calcCutFaces() const
  51 {
  52     if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
  53     {
  54         FatalErrorIn("void enrichedPatch::calcCutFaces() const")
  55             << "Cut faces addressing already calculated."
  56             << abort(FatalError);
  57     }
  58
  59     const faceList& enFaces = enrichedFaces();
  60     const labelList& mp = meshPoints();
  61     const faceList& lf = localFaces();
  62     const pointField& lp = localPoints();
  63     const labelListList& pp = pointPoints();
  64     // Pout<< "enFaces: " << enFaces << endl;
  65     // Pout<< "lf: " << lf << endl;
  66     // Pout<< "lp: " << lp << endl;
  67     // Pout<< "pp: " << pp << endl;
  68     const Map<labelList>& masterPointFaceAddr = masterPointFaces();
  69
  70     // Prepare the storage
  71     DynamicList<face> cf(2*lf.size());
  72     DynamicList<label> cfMaster(2*lf.size());
  73     DynamicList<label> cfSlave(2*lf.size());
  74
  75     // Algorithm
  76     // Go through all the faces
  77     // 1) For each face, start at point zero and grab the edge zero.
  78     //    Both points are added into the cut face.  Mark the first edge
  79     //    as used and position the current point as the end of the first
  80     //    edge and previous point as point zero.
  81     // 2) Grab all possible points from the current point.  Excluding
  82     //    the previous point find the point giving the biggest right
  83     //    turn. Add the point to the face and mark the edges as used.
  84     //    Continue doing this until the face is complete, i.e. the next point
  85     //    to add is the first point of the face.
  86     // 3) Once the facet is completed, register its owner from the face
  87     //    it has been created from (remember that the first lot of faces
  88     //    inserted into the enriched faces list are the slaves, to allow
  89     //    matching of the other side).
  90     // 4) If the facet is created from an original slave face, go
  91     //    through the master patch and try to identify the master face
  92     //    this facet belongs to.  This is recognised by the fact that the
  93     //    faces consists exclusively of the points of the master face and
  94     //    the points projected onto the face.
  95
  96     // Create a set of edge usage parameters
  97     HashSet<edge, Hash<edge> > edgesUsedOnce(pp.size());
  98     HashSet<edge, Hash<edge> > edgesUsedTwice
  99         (pp.size()*primitiveMesh::edgesPerPoint_);
 100
 101
 102     forAll(lf, faceI)
 103     {
 104         const face& curLocalFace = lf[faceI];
 105         const face& curGlobalFace = enFaces[faceI];
 106
 107         // Pout<< "Doing face " << faceI
 108         //     << " local: " << curLocalFace
 109         //     << " or " << curGlobalFace
 110         //     << endl;
 111
 112         // if (faceI < slavePatch_.size())
 113         // {
 114         //     Pout<< "original slave: " << slavePatch_[faceI]
 115         //         << " local: " << slavePatch_.localFaces()[faceI] << endl;
 116         // }
 117         // else
 118         // {
 119         //     Pout<< "original master: "
 120         //         << masterPatch_[faceI - slavePatch_.size()] << " "
 121         //         << masterPatch_.localFaces()[faceI - slavePatch_.size()]
 122         //         << endl;
 123         // }
 124         // {
 125         //     pointField facePoints = curLocalFace.points(lp);
 126         //     forAll(curLocalFace, pointI)
 127         //     {
 128         //         Pout<< "v " << facePoints[pointI].x() << " "
 129         //             << facePoints[pointI].y() << " "
 130         //             << facePoints[pointI].z() << endl;
 131         //     }
 132         // }
 133
 134         // Track the usage of face edges.  When all edges are used, the
 135         // face decomposition is complete.
 136         // The seed edges include all the edges of the original face + all edges
 137         // of other faces that have been used in the construction of the
 138         // facet.  Edges from other faces can be considered as
 139         // internal to the current face if used only once.
 140
 141         // Track the edge usage to avoid duplicate faces and reset it to unused
 142         boolList usedFaceEdges(curLocalFace.size(), false);
 143
 144         SLList<edge> edgeSeeds;
 145
 146         // Insert the edges of current face into the seed list.
 147         edgeList cfe = curLocalFace.edges();
 148         forAll(curLocalFace, edgeI)
 149         {
 150             edgeSeeds.append(cfe[edgeI]);
 151         }
 152
 153         // Grab face normal
 154         vector normal = curLocalFace.normal(lp);
 155         normal /= mag(normal);
 156
 157         while (edgeSeeds.size())
 158         {
 159             // Pout<< "edgeSeeds.size(): "
 160             //     << edgeSeeds.size()
 161             //     << endl;
 162
 163             const edge curEdge = edgeSeeds.removeHead();
 164
 165             // Locate the edge in current face
 166             const label curEdgeWhich = curLocalFace.which(curEdge.start());
 167
 168             // Check if the edge is in current face and if it has been
 169             // used already.  If so, skip it
 170             if
 171             (
 172                 curEdgeWhich > -1
 173              && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
 174             )
 175             {
 176                 // Edge is in the starting face.
 177                 // If unused, mark it as used; if used, skip it
 178                 if (usedFaceEdges[curEdgeWhich]) continue;
 179
 180                 usedFaceEdges[curEdgeWhich] = true;
 181             }
 182
 183             // If the edge has already been used twice, skip it
 184             if (edgesUsedTwice.found(curEdge)) continue;
 185
 186             // Pout<< "Trying new edge (" << mp[curEdge.start()]
 187             //     << ", " << mp[curEdge.end()]
 188             //     << ") seed: " << curEdge
 189             //     << " used: " << edgesUsedTwice.found(curEdge)
 190             //     << endl;
 191
 192             // Estimate the size of cut face as twice the size of original face
 193             DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
 194             DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
 195
 196             // Found unused edge.
 197             label prevPointLabel = curEdge.start();
 198             cutFaceGlobalPoints.append(mp[prevPointLabel]);
 199             cutFaceLocalPoints.append(prevPointLabel);
 200             // Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;
 201             // Grab current point and append it to the list
 202             label curPointLabel = curEdge.end();
 203             point curPoint = lp[curPointLabel];
 204             cutFaceGlobalPoints.append(mp[curPointLabel]);
 205             cutFaceLocalPoints.append(curPointLabel);
 206
 207             bool completedCutFace = false;
 208
 209             label faceSizeDebug = maxFaceSizeDebug_;
 210
 211             do
 212             {
 213                 // Grab the next point options
 214
 215                 // Pout<< "curPointLabel: " << mp[curPointLabel] << endl;
 216
 217                 const labelList& nextPoints = pp[curPointLabel];
 218
 219                 // Pout<< "nextPoints: "
 220                 //     << UIndirectList<label>(mp, nextPoints)
 221                 //     << endl;
 222
 223                 // Get the vector along the edge and the right vector
 224                 vector ahead = curPoint - lp[prevPointLabel];
 225                 ahead -= normal*(normal & ahead);
 226                 ahead /= mag(ahead);
 227
 228                 vector right = normal ^ ahead;
 229                 right /= mag(right);
 230
 231                 // Pout<< "normal: " << normal
 232                 //     << " ahead: " << ahead
 233                 //     << " right: " << right
 234                 //     << endl;
 235
 236                 scalar atanTurn = -GREAT;
 237                 label bestAtanPoint = -1;
 238
 239                 forAll(nextPoints, nextI)
 240                 {
 241                     // Exclude the point we are coming from; there will always
 242                     // be more than one edge, so this is safe
 243                     if (nextPoints[nextI] != prevPointLabel)
 244                     {
 245                         // Pout<< "cur point: " << curPoint
 246                         //     << " trying for point: "
 247                         //     << mp[nextPoints[nextI]]
 248                         //     << " " << lp[nextPoints[nextI]];
 249                         vector newDir = lp[nextPoints[nextI]] - curPoint;
 250                         // Pout<< " newDir: " << newDir
 251                         //     << " mag: " << mag(newDir) << flush;
 252                         newDir -= normal*(normal & newDir);
 253                         scalar magNewDir = mag(newDir);
 254                         // Pout<< " corrected: " << newDir
 255                         //     << " mag: " << mag(newDir) << flush;
 256
 257                         if (magNewDir < SMALL)
 258                         {
 259                             FatalErrorIn
 260                             (
 261                                 "void enrichedPatch::"
 262                                 "calcCutFaces() const"
 263                             )   << "Zero length edge detected.  Probable "
 264                                 << "projection error: slave patch probably "
 265                                 << "does not project onto master.  "
 266                                 << "Please switch on "
 267                                 << "enriched patch debug for more info"
 268                                 << abort(FatalError);
 269                         }
 270
 271                         newDir /= magNewDir;
 272
 273                         scalar curAtanTurn =
 274                             atan2(newDir & right, newDir & ahead);
 275
 276                         // Pout<< " atan: " << curAtanTurn << endl;
 277
 278                         if (curAtanTurn > atanTurn)
 279                         {
 280                             bestAtanPoint = nextPoints[nextI];
 281                             atanTurn = curAtanTurn;
 282                         }
 283                     } // end of prev point skip
 284                 } // end of next point selection
 285
 286                 // Pout<< "   bestAtanPoint: " << bestAtanPoint << " or "
 287                 //     << mp[bestAtanPoint]
 288                 //     << endl;
 289
 290                 // Selected next best point.
 291
 292                 // Pout<< "cutFaceGlobalPoints: "
 293                 //     << cutFaceGlobalPoints
 294                 //     << endl;
 295
 296                 // Check if the edge about to be added has been used
 297                 // in the current face or twice in other faces.  If
 298                 // so, the face is bad.
 299                 const label whichNextPoint = curLocalFace.which(curPointLabel);
 300
 301                 if
 302                 (
 303                     whichNextPoint > -1
 304                  && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
 305                  && usedFaceEdges[whichNextPoint]
 306                 )
 307                 {
 308                     // This edge is already used in current face
 309                     // face cannot be good; start on a new one
 310
 311                     // Pout<< "Double usage in current face, cannot be good"
 312                     //     << endl;
 313
 314                     completedCutFace = true;
 315                 }
 316
 317
 318                 if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
 319                 {
 320                     // This edge is already used -
 321                     // face cannot be good; start on a new one
 322                     completedCutFace = true;
 323
 324                     // Pout<< "Double usage elsewhere, cannot be good" << endl;
 325                 }
 326
 327                 if (completedCutFace) continue;
 328
 329                 // Insert the next best point into the list
 330                 if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
 331                 {
 332                     // Face is completed.  Save it and move on to the next
 333                     // available edge
 334                     completedCutFace = true;
 335
 336                     if (debug)
 337                     {
 338                         Pout<< " local: " << cutFaceLocalPoints
 339                             << " one side: " << faceI;
 340                     }
 341
 342                     // Append the face
 343                     face cutFaceGlobal;
 344                     cutFaceGlobal.transfer(cutFaceGlobalPoints);
 345
 346                     cf.append(cutFaceGlobal);
 347
 348                     face cutFaceLocal;
 349                     cutFaceLocal.transfer(cutFaceLocalPoints);
 350
 351                     // Pout<< "\ncutFaceLocal: "
 352                     //     << cutFaceLocal.points(lp)
 353                     //     << endl;
 354
 355                     // Go through all edges of the cut faces.
 356                     // If the edge corresponds to a starting face edge,
 357                     // mark the starting face edge as true
 358
 359                     forAll(cutFaceLocal, cutI)
 360                     {
 361                         const edge curCutFaceEdge
 362                         (
 363                             cutFaceLocal[cutI],
 364                             cutFaceLocal.nextLabel(cutI)
 365                         );
 366
 367                         // Increment the usage count using two hash sets
 368                         HashSet<edge, Hash<edge> >::iterator euoIter =
 369                             edgesUsedOnce.find(curCutFaceEdge);
 370
 371                         if (euoIter == edgesUsedOnce.end())
 372                         {
 373                             // Pout<< "Found edge not used before: "
 374                             //     << curCutFaceEdge
 375                             //     << endl;
 376                             edgesUsedOnce.insert(curCutFaceEdge);
 377                         }
 378                         else
 379                         {
 380                             // Pout<< "Found edge used once: "
 381                             //     << curCutFaceEdge
 382                             //     << endl;
 383                             edgesUsedOnce.erase(euoIter);
 384                             edgesUsedTwice.insert(curCutFaceEdge);
 385                         }
 386
 387                         const label curCutFaceEdgeWhich = curLocalFace.which
 388                         (
 389                             curCutFaceEdge.start()
 390                         );
 391
 392                         if
 393                         (
 394                             curCutFaceEdgeWhich > -1
 395                          && curLocalFace.nextLabel(curCutFaceEdgeWhich)
 396                          == curCutFaceEdge.end()
 397                         )
 398                         {
 399                             // Found edge in original face
 400
 401                             // Pout<< "Found edge in orig face: "
 402                             //     << curCutFaceEdge << ": "
 403                             //     << curCutFaceEdgeWhich
 404                             //     << endl;
 405
 406                             usedFaceEdges[curCutFaceEdgeWhich] = true;
 407                         }
 408                         else
 409                         {
 410                             // Edge not in original face.  Add it to seeds
 411
 412                             // Pout<< "Found new edge seed: "
 413                             //     << curCutFaceEdge
 414                             //     << endl;
 415
 416                             edgeSeeds.append(curCutFaceEdge.reverseEdge());
 417                         }
 418                     }
 419
 420                     // Find out what the other side is
 421
 422                     // Algorithm
 423
 424                     // If the face is in the slave half of the
 425                     // enrichedFaces lists, it may be matched against
 426                     // the master face.  It can be recognised by the
 427                     // fact that all its points belong to one master
 428                     // face.  For this purpose:
 429                     // 1) Grab the master faces around the point zero
 430                     // of the cut face and collect all master faces
 431                     // around it.
 432                     // 2) Loop through the rest of cut face points and
 433                     // if the point refers to one of the faces marked
 434                     // by point zero, increment its count.
 435                     // 3) When completed, the face whose count is
 436                     // equal to the number of points in the cut face
 437                     // is the other side.  If this is not the case, there is no
 438                     // face on the other side.
 439
 440                     if (faceI < slavePatch_.size())
 441                     {
 442                         Map<labelList>::const_iterator mpfAddrIter =
 443                             masterPointFaceAddr.find(cutFaceGlobal[0]);
 444
 445                         bool otherSideFound = false;
 446
 447                         if (mpfAddrIter != masterPointFaceAddr.end())
 448                         {
 449                             bool miss = false;
 450
 451                             // Create the label count using point zero
 452                             const labelList& masterFacesOfPZero = mpfAddrIter();
 453
 454                             labelList hits(masterFacesOfPZero.size(), 1);
 455
 456                             for
 457                             (
 458                                 label pointI = 1;
 459                                 pointI < cutFaceGlobal.size();
 460                                 pointI++
 461                             )
 462                             {
 463                                 Map<labelList>::const_iterator
 464                                     mpfAddrPointIter =
 465                                         masterPointFaceAddr.find
 466                                         (
 467                                             cutFaceGlobal[pointI]
 468                                         );
 469
 470                                 if
 471                                 (
 472                                     mpfAddrPointIter
 473                                  == masterPointFaceAddr.end()
 474                                 )
 475                                 {
 476                                     // Point is off the master patch. Skip
 477                                     miss = true;
 478                                     break;
 479                                 }
 480
 481                                 const labelList& curMasterFaces =
 482                                     mpfAddrPointIter();
 483
 484                                 // For every current face, try to find it in the
 485                                 // zero-list
 486                                 forAll(curMasterFaces, i)
 487                                 {
 488                                     forAll(masterFacesOfPZero, j)
 489                                     {
 490                                         if
 491                                         (
 492                                             curMasterFaces[i]
 493                                          == masterFacesOfPZero[j]
 494                                         )
 495                                         {
 496                                             hits[j]++;
 497                                             break;
 498                                         }
 499                                     }
 500                                 }
 501                             }
 502
 503                             // If all point are found attempt matching
 504                             if (!miss)
 505                             {
 506                                 forAll(hits, pointI)
 507                                 {
 508                                     if (hits[pointI] == cutFaceGlobal.size())
 509                                     {
 510                                         // Found other side.
 511                                         otherSideFound = true;
 512
 513                                         cfMaster.append
 514                                         (
 515                                             masterFacesOfPZero[pointI]
 516                                         );
 517
 518                                         cfSlave.append(faceI);
 519
 520                                         // Reverse the face such that it
 521                                         // points out of the master patch
 522                                         cf.last().flip();
 523
 524                                         if (debug)
 525                                         {
 526                                             Pout<< " other side: "
 527                                                 << masterFacesOfPZero[pointI]
 528                                                 << endl;
 529                                         }
 530                                     } // end of hits
 531                                 } // end of for all hits
 532
 533                             } // end of not miss
 534
 535                             if (!otherSideFound || miss)
 536                             {
 537                                 if (debug)
 538                                 {
 539                                     Pout<< " solo slave A" << endl;
 540                                 }
 541
 542                                 cfMaster.append(-1);
 543                                 cfSlave.append(faceI);
 544                             }
 545                         }
 546                         else
 547                         {
 548                             // First point not in master patch
 549                             if (debug)
 550                             {
 551                                 Pout<< " solo slave B" << endl;
 552                             }
 553
 554                             cfMaster.append(-1);
 555                             cfSlave.append(faceI);
 556                         }
 557                     }
 558                     else
 559                     {
 560                         if (debug)
 561                         {
 562                             Pout<< " master side" << endl;
 563                         }
 564
 565                         cfMaster.append(faceI - slavePatch_.size());
 566                         cfSlave.append(-1);
 567                     }
 568                 }
 569                 else
 570                 {
 571                     // Face not completed.  Prepare for the next point search
 572                     prevPointLabel = curPointLabel;
 573                     curPointLabel = bestAtanPoint;
 574                     curPoint = lp[curPointLabel];
 575                     cutFaceGlobalPoints.append(mp[curPointLabel]);
 576                     cutFaceLocalPoints.append(curPointLabel);
 577
 578                     if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
 579                     {
 580                         faceSizeDebug *= 2;
 581
 582                         // Check for duplicate points in the face
 583                         forAll(cutFaceGlobalPoints, checkI)
 584                         {
 585                             for
 586                             (
 587                                 label checkJ = checkI + 1;
 588                                 checkJ < cutFaceGlobalPoints.size();
 589                                 checkJ++
 590                             )
 591                             {
 592                                 if
 593                                 (
 594                                     cutFaceGlobalPoints[checkI]
 595                                  == cutFaceGlobalPoints[checkJ]
 596                                 )
 597                                 {
 598                                     // Shrink local points for debugging
 599                                     cutFaceLocalPoints.shrink();
 600
 601                                     face origFace;
 602                                     face origFaceLocal;
 603                                     if (faceI < slavePatch_.size())
 604                                     {
 605                                         origFace = slavePatch_[faceI];
 606                                         origFaceLocal =
 607                                             slavePatch_.localFaces()[faceI];
 608                                     }
 609                                     else
 610                                     {
 611                                         origFace =
 612                                             masterPatch_
 613                                             [faceI - slavePatch_.size()];
 614
 615                                         origFaceLocal =
 616                                             masterPatch_.localFaces()
 617                                             [faceI - slavePatch_.size()];
 618                                     }
 619
 620                                     FatalErrorIn
 621                                     (
 622                                         "void enrichedPatch::"
 623                                         "calcCutFaces() const"
 624                                     )   << "Duplicate point found in cut face. "
 625                                         << "Error in the face cutting "
 626                                         << "algorithm for global face "
 627                                         << origFace << " local face "
 628                                         << origFaceLocal << nl
 629                                         << "Slave size: " << slavePatch_.size()
 630                                         << " Master size: "
 631                                         << masterPatch_.size()
 632                                         << " index: " << faceI << ".\n"
 633                                         << "Face: " << curGlobalFace << nl
 634                                         << "Cut face: " << cutFaceGlobalPoints
 635                                         << " local: " << cutFaceLocalPoints
 636                                         << nl << "Points: "
 637                                         << face(cutFaceLocalPoints).points(lp)
 638                                         << abort(FatalError);
 639                                 }
 640                             }
 641                         }
 642                     } // end of debug
 643                 }
 644             } while (!completedCutFace);
 645         } // end of used edges
 646
 647         if (debug)
 648         {
 649             Pout<< " Finished face " << faceI << endl;
 650         }
 651
 652     } // end of local faces
 653
 654     // Re-pack the list into compact storage
 655     cutFacesPtr_ = new faceList();
 656     cutFacesPtr_->transfer(cf);
 657
 658     cutFaceMasterPtr_ = new labelList();
 659     cutFaceMasterPtr_->transfer(cfMaster);
 660
 661     cutFaceSlavePtr_ = new labelList();
 662     cutFaceSlavePtr_->transfer(cfSlave);
 663 }
 664
 665
 666 void Foam::enrichedPatch::clearCutFaces()
 667 {
 668     deleteDemandDrivenData(cutFacesPtr_);
 669     deleteDemandDrivenData(cutFaceMasterPtr_);
 670     deleteDemandDrivenData(cutFaceSlavePtr_);
 671 }
 672
 673
 674 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 675
 676 const Foam::faceList& Foam::enrichedPatch::cutFaces() const
 677 {
 678     if (!cutFacesPtr_)
 679     {
 680         calcCutFaces();
 681     }
 682
 683     return *cutFacesPtr_;
 684 }
 685
 686
 687 const Foam::labelList& Foam::enrichedPatch::cutFaceMaster() const
 688 {
 689     if (!cutFaceMasterPtr_)
 690     {
 691         calcCutFaces();
 692     }
 693
 694     return *cutFaceMasterPtr_;
 695 }
 696
 697
 698 const Foam::labelList& Foam::enrichedPatch::cutFaceSlave() const
 699 {
 700     if (!cutFaceSlavePtr_)
 701     {
 702         calcCutFaces();
 703     }
 704
 705     return *cutFaceSlavePtr_;
 706 }
 707
 708
 709 // ************************************************************************* //