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