src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.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 "slidingInterface.H"
  28 #include "polyMesh.H"
  29 #include "primitiveMesh.H"
  30 #include "line.H"
  31 #include "triPointRef.H"
  32 #include "plane.H"
  33 #include "polyTopoChanger.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 const Foam::scalar Foam::slidingInterface::pointMergeTol_ = 0.05;
  38 const Foam::scalar Foam::slidingInterface::edgeMergeTol_ = 0.01;
  39 const Foam::label Foam::slidingInterface::nFacesPerSlaveEdge_ = 5;
  40 const Foam::label Foam::slidingInterface::edgeFaceEscapeLimit_ = 10;
  41
  42 const Foam::scalar Foam::slidingInterface::integralAdjTol_ = 0.05;
  43 const Foam::scalar Foam::slidingInterface::edgeMasterCatchFraction_ = 0.4;
  44 const Foam::scalar Foam::slidingInterface::edgeEndCutoffTol_ = 0.0001;
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Index of debug signs:
  50 // a - integral match adjustment: point adjusted
  51 // n - integral match adjustment: point not adjusted
  52 // . - loop of the edge-to-face interaction detection
  53 // x - reversal of direction in edge-to-face interaction detection
  54 // + - complete edge-to-face interaction detection
  55 // z - incomplete edge-to-face interaction detection.  This may be OK for edges
  56 //     crossing from one to the other side of multiply connected master patch
  57 // * - colinear triangle: adjusting projection with slave face normal
  58 // m - master point inserted into the edge
  59
  60 bool Foam::slidingInterface::projectPoints() const
  61 {
  62     if (debug)
  63     {
  64         Pout<< "bool slidingInterface::projectPoints() : "
  65             << " for object " << name() << " : "
  66             << "Projecting slave points onto master surface." << endl;
  67     }
  68
  69     // Algorithm:
  70     // 1) Go through all the points of the master and slave patch and calculate
  71     //    minimum edge length coming from the point.  Calculate the point
  72     //    merge tolerance as the fraction of mimimum edge length.
  73     // 2) Project all the slave points onto the master patch
  74     //    in the normal direction.
  75     // 3) If some points have missed and the match is integral, the
  76     //    points in question will be adjusted.  Find the nearest face for
  77     //    those points and continue with the procedure.
  78     // 4) For every point, check all the points of the face it has hit.
  79     //    For every pair of points find if their distance is smaller than
  80     //    both the master and slave merge tolerance.  If so, the slave point
  81     //    is moved to the location of the master point.  Remember the master
  82     //    point index.
  83     // 5) For every unmerged slave point, check its distance to all the
  84     //    edges of the face it has hit.  If the distance is smaller than the
  85     //    edge merge tolerance, the point will be moved onto the edge.
  86     //    Remember the master edge index.
  87     // 6) The remaning slave points will be projected into faces.  Remember the
  88     //    master face index.
  89     // 7) For the points that miss the master patch, grab the nearest face
  90     //    on the master and leave the slave point where it started
  91     //    from and the miss is recorded.
  92
  93     const polyMesh& mesh = topoChanger().mesh();
  94
  95     const primitiveFacePatch& masterPatch =
  96         mesh.faceZones()[masterFaceZoneID_.index()]();
  97
  98     const primitiveFacePatch& slavePatch =
  99         mesh.faceZones()[slaveFaceZoneID_.index()]();
 100
 101     // Get references to local points, local edges and local faces
 102     // for master and slave patch
 103 //     const labelList& masterMeshPoints = masterPatch.meshPoints();
 104     const pointField& masterLocalPoints = masterPatch.localPoints();
 105     const faceList& masterLocalFaces = masterPatch.localFaces();
 106     const edgeList& masterEdges = masterPatch.edges();
 107     const labelListList& masterEdgeFaces = masterPatch.edgeFaces();
 108     const labelListList& masterFaceEdges = masterPatch.faceEdges();
 109     const labelListList& masterFaceFaces = masterPatch.faceFaces();
 110 //     Pout<< "Master patch.  Local points: " << masterLocalPoints << nl
 111 //         << "Master patch.  Mesh points: " << masterPatch.meshPoints() << nl
 112 //         << "Local faces: " << masterLocalFaces << nl
 113 //         << "local edges: " << masterEdges << endl;
 114
 115 //     const labelList& slaveMeshPoints = slavePatch.meshPoints();
 116     const pointField& slaveLocalPoints = slavePatch.localPoints();
 117     const edgeList& slaveEdges = slavePatch.edges();
 118     const labelListList& slaveEdgeFaces = slavePatch.edgeFaces();
 119     const vectorField& slavePointNormals = slavePatch.pointNormals();
 120 //     const vectorField& slaveFaceNormals = slavePatch.faceNormals();
 121 //     Pout<< "Slave patch.  Local points: " << slaveLocalPoints << nl
 122 //         << "Slave patch.  Mesh points: " << slavePatch.meshPoints() << nl
 123 //         << "Local faces: " << slavePatch.localFaces() << nl
 124 //         << "local edges: " << slaveEdges << endl;
 125
 126     // Calculate min edge distance for points and faces
 127
 128     // Calculate min edge length for the points and faces of master patch
 129     scalarField minMasterPointLength(masterLocalPoints.size(), GREAT);
 130     scalarField minMasterFaceLength(masterPatch.size(), GREAT);
 131
 132     forAll (masterEdges, edgeI)
 133     {
 134         const edge& curEdge = masterEdges[edgeI];
 135
 136         const scalar curLength =
 137             masterEdges[edgeI].mag(masterLocalPoints);
 138
 139         // Do points
 140         minMasterPointLength[curEdge.start()] =
 141             min
 142             (
 143                 minMasterPointLength[curEdge.start()],
 144                 curLength
 145             );
 146
 147         minMasterPointLength[curEdge.end()] =
 148             min
 149             (
 150                 minMasterPointLength[curEdge.end()],
 151                 curLength
 152             );
 153
 154         // Do faces
 155         const labelList& curFaces = masterEdgeFaces[edgeI];
 156
 157         forAll (curFaces, faceI)
 158         {
 159             minMasterFaceLength[curFaces[faceI]] =
 160                 min
 161                 (
 162                     minMasterFaceLength[curFaces[faceI]],
 163                     curLength
 164                 );
 165         }
 166     }
 167
 168 //     Pout << "min length for master points: " << minMasterPointLength << endl
 169 //         << "min length for master faces: " << minMasterFaceLength << endl;
 170
 171     // Calculate min edge length for the points and faces of slave patch
 172     scalarField minSlavePointLength(slaveLocalPoints.size(), GREAT);
 173     scalarField minSlaveFaceLength(slavePatch.size(), GREAT);
 174
 175     forAll (slaveEdges, edgeI)
 176     {
 177         const edge& curEdge = slaveEdges[edgeI];
 178
 179         const scalar curLength =
 180             slaveEdges[edgeI].mag(slaveLocalPoints);
 181
 182         // Do points
 183         minSlavePointLength[curEdge.start()] =
 184             min
 185             (
 186                 minSlavePointLength[curEdge.start()],
 187                 curLength
 188             );
 189
 190         minSlavePointLength[curEdge.end()] =
 191             min
 192             (
 193                 minSlavePointLength[curEdge.end()],
 194                 curLength
 195             );
 196
 197         // Do faces
 198         const labelList& curFaces = slaveEdgeFaces[edgeI];
 199
 200         forAll (curFaces, faceI)
 201         {
 202             minSlaveFaceLength[curFaces[faceI]] =
 203                 min
 204                 (
 205                     minSlaveFaceLength[curFaces[faceI]],
 206                     curLength
 207                 );
 208         }
 209     }
 210
 211 //     Pout << "min length for slave points: " << minSlavePointLength << endl
 212 //         << "min length for slave faces: " << minSlaveFaceLength << endl;
 213
 214     // Project slave points onto the master patch
 215
 216     // Face hit by the slave point
 217     List<objectHit> slavePointFaceHits =
 218         slavePatch.projectPoints
 219         (
 220             masterPatch,
 221             slavePointNormals,
 222             projectionAlgo_
 223         );
 224
 225 //     Pout << "USING N-SQAURED!!!" << endl;
 226 //     List<objectHit> slavePointFaceHits =
 227 //         projectPointsNSquared<face, List, const pointField&>
 228 //         (
 229 //             slavePatch,
 230 //             masterPatch,
 231 //             slavePointNormals,
 232 //             projectionAlgo_
 233 //         );
 234
 235     if (debug)
 236     {
 237         label nHits = 0;
 238
 239         forAll (slavePointFaceHits, pointI)
 240         {
 241             if (slavePointFaceHits[pointI].hit())
 242             {
 243                 nHits++;
 244             }
 245         }
 246
 247         Pout<< "Number of hits in point projection: " << nHits
 248             << " out of " << slavePointFaceHits.size() << " points."
 249             << endl;
 250     }
 251
 252     // Projected slave points are stored for mesh motion correction
 253     if (projectedSlavePointsPtr_) delete projectedSlavePointsPtr_;
 254
 255     projectedSlavePointsPtr_ =
 256         new pointField(slavePointFaceHits.size(), vector::zero);
 257     pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
 258
 259     // Adjust projection to type of match
 260
 261     label nAdjustedPoints = 0;
 262
 263     // If the match is integral and some points have missed,
 264     // find nearest master face
 265     if (matchType_ == INTEGRAL)
 266     {
 267         if (debug)
 268         {
 269             Pout<< "bool slidingInterface::projectPoints() for object "
 270                 << name() << " : "
 271                 << "Adjusting point projection for integral match: ";
 272         }
 273
 274         forAll (slavePointFaceHits, pointI)
 275         {
 276             if (slavePointFaceHits[pointI].hit())
 277             {
 278                 // Grab the hit point
 279                 projectedSlavePoints[pointI] =
 280                     masterLocalFaces
 281                         [slavePointFaceHits[pointI].hitObject()].ray
 282                         (
 283                             slaveLocalPoints[pointI],
 284                             slavePointNormals[pointI],
 285                             masterLocalPoints,
 286                             projectionAlgo_
 287                         ).hitPoint();
 288             }
 289             else
 290             {
 291                 // Grab the nearest point on the face (edge)
 292                 pointHit missAdjust =
 293                     masterLocalFaces[slavePointFaceHits[pointI].hitObject()].ray
 294                     (
 295                         slaveLocalPoints[pointI],
 296                         slavePointNormals[pointI],
 297                         masterLocalPoints,
 298                         projectionAlgo_
 299                     );
 300
 301                 const point nearPoint = missAdjust.missPoint();
 302                 const point missPoint =
 303                     slaveLocalPoints[pointI]
 304                   + slavePointNormals[pointI]*missAdjust.distance();
 305
 306                 // Calculate the tolerance
 307                 const scalar mergeTol =
 308                     integralAdjTol_*minSlavePointLength[pointI];
 309
 310                 // Adjust the hit
 311                 if (mag(nearPoint - missPoint) < mergeTol)
 312                 {
 313                     if (debug)
 314                     {
 315                         Pout << "a";
 316                     }
 317
 318 //                     Pout<< "Moving slave point in integral adjustment "
 319 //                         << pointI << " miss point: " << missPoint
 320 //                         << " near point: " << nearPoint
 321 //                         << " mergeTol: " << mergeTol
 322 //                         << " dist: " << mag(nearPoint - missPoint) << endl;
 323
 324                     projectedSlavePoints[pointI] = nearPoint;
 325
 326                     slavePointFaceHits[pointI] =
 327                         objectHit(true, slavePointFaceHits[pointI].hitObject());
 328
 329                     nAdjustedPoints++;
 330                 }
 331                 else
 332                 {
 333                     projectedSlavePoints[pointI] = slaveLocalPoints[pointI];
 334
 335                     if (debug)
 336                     {
 337                         Pout << "n";
 338                     }
 339                 }
 340             }
 341         }
 342
 343         if (debug)
 344         {
 345             Pout << " done." << endl;
 346         }
 347     }
 348     else if (matchType_ == PARTIAL)
 349     {
 350         forAll (slavePointFaceHits, pointI)
 351         {
 352             if (slavePointFaceHits[pointI].hit())
 353             {
 354                 // Grab the hit point
 355                 projectedSlavePoints[pointI] =
 356                     masterLocalFaces
 357                         [slavePointFaceHits[pointI].hitObject()].ray
 358                         (
 359                             slaveLocalPoints[pointI],
 360                             slavePointNormals[pointI],
 361                             masterLocalPoints,
 362                             projectionAlgo_
 363                         ).hitPoint();
 364             }
 365             else
 366             {
 367                 // The point remains where it started from
 368                 projectedSlavePoints[pointI] = slaveLocalPoints[pointI];
 369             }
 370         }
 371     }
 372     else
 373     {
 374         FatalErrorIn
 375         (
 376             "bool slidingInterface::projectPoints() const"
 377         )   << "Problem in point projection.  Unknown sliding match type "
 378             << " for object " << name()
 379             << abort(FatalError);
 380     }
 381
 382     if (debug)
 383     {
 384         Pout<< "Number of adjusted points in projection: "
 385             << nAdjustedPoints << endl;
 386
 387         // Check for zero-length edges in slave projection
 388         scalar minEdgeLength = GREAT;
 389         scalar el = 0;
 390         label nShortEdges = 0;
 391
 392         forAll (slaveEdges, edgeI)
 393         {
 394             el = slaveEdges[edgeI].mag(projectedSlavePoints);
 395
 396             if (el < SMALL)
 397             {
 398                 Pout<< "Point projection problems for edge: "
 399                     << slaveEdges[edgeI] << ". Length = " << el
 400                     << endl;
 401
 402                 nShortEdges++;
 403             }
 404
 405             minEdgeLength = min(minEdgeLength, el);
 406         }
 407
 408         if (nShortEdges > 0)
 409         {
 410             FatalErrorIn
 411             (
 412                 "bool slidingInterface::projectPoints() const"
 413             )   << "Problem in point projection.  " << nShortEdges
 414                 << " short projected edges "
 415                 << "after adjustment for object " << name()
 416                 << abort(FatalError);
 417         }
 418         else
 419         {
 420             Pout << " ... projection OK." << endl;
 421         }
 422     }
 423 //     scalarField magDiffs(mag(slaveLocalPoints - projectedSlavePoints));
 424 //     Pout<< "Slave point face hits: " << slavePointFaceHits << nl
 425 //         << "slave points: " << slaveLocalPoints << nl
 426 //         << "Projected slave points: " << projectedSlavePoints
 427 //         << "diffs: " << magDiffs << endl;
 428
 429     // Merge projected points against master points
 430
 431     labelList slavePointPointHits(slaveLocalPoints.size(), -1);
 432     labelList masterPointPointHits(masterLocalPoints.size(), -1);
 433
 434     // Go through all the slave points and compare them against all the points
 435     // in the master face they hit.  If the distance between the projected point
 436     // and any of the master face points is smaller than both tolerances,
 437     // merge the projected point by:
 438     // 1) adjusting the projected point to coincide with the
 439     //    master point it merges with
 440     // 2) remembering the hit master point index in slavePointPointHits
 441     // 3) resetting the hit face to -1
 442     // 4) If a master point has been hit directly, it cannot be merged
 443     // into the edge. Mark it as used in the reverse list
 444
 445     label nMergedPoints = 0;
 446
 447     forAll (projectedSlavePoints, pointI)
 448     {
 449         if (slavePointFaceHits[pointI].hit())
 450         {
 451             // Taking a non-const reference so the point can be adjusted
 452             point& curPoint = projectedSlavePoints[pointI];
 453
 454             // Get the hit face
 455             const face& hitFace =
 456                 masterLocalFaces[slavePointFaceHits[pointI].hitObject()];
 457
 458             label mergePoint = -1;
 459             scalar mergeDist = GREAT;
 460
 461             // Try all point before deciding on best fit.
 462             forAll (hitFace, hitPointI)
 463             {
 464                 scalar dist =
 465                     mag(masterLocalPoints[hitFace[hitPointI]] - curPoint);
 466
 467                 // Calculate the tolerance
 468                 const scalar mergeTol =
 469                     pointMergeTol_*
 470                     min
 471                     (
 472                         minSlavePointLength[pointI],
 473                         minMasterPointLength[hitFace[hitPointI]]
 474                     );
 475
 476                 if (dist < mergeTol && dist < mergeDist)
 477                 {
 478                     mergePoint = hitFace[hitPointI];
 479                     mergeDist = dist;
 480
 481 //                     Pout<< "Merging slave point "
 482 //                         << slavePatch.meshPoints()[pointI] << " at "
 483 //                         << slaveLocalPoints[pointI] << " with master "
 484 //                         << masterPatch.meshPoints()[mergePoint] << " at "
 485 //                         << masterLocalPoints[mergePoint]
 486 //                         << ". dist: " << mergeDist
 487 //                         << " mergeTol: " << mergeTol << endl;
 488                 }
 489             }
 490
 491             if (mergePoint > -1)
 492             {
 493                 // Point is to be merged with master point
 494                 nMergedPoints++;
 495
 496                 slavePointPointHits[pointI] = mergePoint;
 497                 curPoint = masterLocalPoints[mergePoint];
 498                 masterPointPointHits[mergePoint] = pointI;
 499             }
 500         }
 501     }
 502
 503 //     Pout<< "slavePointPointHits: " << slavePointPointHits << nl
 504 //         << "masterPointPointHits: " << masterPointPointHits << endl;
 505
 506     if (debug)
 507     {
 508         // Check for zero-length edges in slave projection
 509         scalar minEdgeLength = GREAT;
 510         scalar el = 0;
 511
 512         forAll (slaveEdges, edgeI)
 513         {
 514             el = slaveEdges[edgeI].mag(projectedSlavePoints);
 515
 516             if (el < SMALL)
 517             {
 518                 Pout<< "Point projection problems for edge: "
 519                     << slaveEdges[edgeI] << ". Length = " << el
 520                     << endl;
 521             }
 522
 523             minEdgeLength = min(minEdgeLength, el);
 524         }
 525
 526         if (minEdgeLength < SMALL)
 527         {
 528             FatalErrorIn
 529             (
 530                 "bool slidingInterface::projectPoints() const"
 531             )   << "Problem in point projection.  Short projected edge "
 532                 << " after point merge for object " << name()
 533                 << abort(FatalError);
 534         }
 535         else
 536         {
 537             Pout << " ... point merge OK." << endl;
 538         }
 539     }
 540
 541     // Merge projected points against master edges
 542
 543     labelList slavePointEdgeHits(slaveLocalPoints.size(), -1);
 544
 545     label nMovedPoints = 0;
 546
 547     forAll (projectedSlavePoints, pointI)
 548     {
 549         // Eliminate the points merged into points
 550         if (slavePointPointHits[pointI] < 0)
 551         {
 552             // Get current point position
 553             point& curPoint = projectedSlavePoints[pointI];
 554
 555             // Get the hit face
 556             const labelList& hitFaceEdges =
 557                 masterFaceEdges[slavePointFaceHits[pointI].hitObject()];
 558
 559             // Calculate the tolerance
 560             const scalar mergeLength =
 561                 min
 562                 (
 563                     minSlavePointLength[pointI],
 564                     minMasterFaceLength[slavePointFaceHits[pointI].hitObject()]
 565                 );
 566
 567             const scalar mergeTol = pointMergeTol_*mergeLength;
 568
 569             scalar minDistance = GREAT;
 570
 571             forAll (hitFaceEdges, edgeI)
 572             {
 573                 const edge& curEdge = masterEdges[hitFaceEdges[edgeI]];
 574
 575                 pointHit edgeHit =
 576                     curEdge.line(masterLocalPoints).nearestDist(curPoint);
 577
 578                 if (edgeHit.hit())
 579                 {
 580                     scalar dist =
 581                         mag(edgeHit.hitPoint() - projectedSlavePoints[pointI]);
 582
 583                     if (dist < mergeTol && dist < minDistance)
 584                     {
 585                         // Point is to be moved onto master edge
 586                         nMovedPoints++;
 587
 588                         slavePointEdgeHits[pointI] = hitFaceEdges[edgeI];
 589                         projectedSlavePoints[pointI] = edgeHit.hitPoint();
 590
 591                         minDistance = dist;
 592
 593 //                         Pout<< "Moving slave point "
 594 //                             << slavePatch.meshPoints()[pointI]
 595 //                             << " (" << pointI
 596 //                             << ") at " << slaveLocalPoints[pointI]
 597 //                             << " onto master edge " << hitFaceEdges[edgeI]
 598 //                             << " or ("
 599 //                             << masterLocalPoints[curEdge.start()]
 600 //                             << masterLocalPoints[curEdge.end()]
 601 //                             << ") hit: " << edgeHit.hitPoint()
 602 //                             << ". dist: " << dist
 603 //                             << " mergeTol: " << mergeTol << endl;
 604                     }
 605                 }
 606             } // end of hit face edges
 607
 608             if (slavePointEdgeHits[pointI] > -1)
 609             {
 610                 // Projected slave point has moved.  Re-attempt merge with
 611                 // master points of the edge
 612                 point& curPoint = projectedSlavePoints[pointI];
 613
 614                 const edge& hitMasterEdge =
 615                     masterEdges[slavePointEdgeHits[pointI]];
 616
 617                 label mergePoint = -1;
 618                 scalar mergeDist = GREAT;
 619
 620                 forAll (hitMasterEdge, hmeI)
 621                 {
 622                     scalar hmeDist =
 623                         mag(masterLocalPoints[hitMasterEdge[hmeI]] - curPoint);
 624
 625                     // Calculate the tolerance
 626                     const scalar mergeTol =
 627                         pointMergeTol_*
 628                         min
 629                         (
 630                             minSlavePointLength[pointI],
 631                             minMasterPointLength[hitMasterEdge[hmeI]]
 632                     );
 633
 634                     if (hmeDist < mergeTol && hmeDist < mergeDist)
 635                     {
 636                         mergePoint = hitMasterEdge[hmeI];
 637                         mergeDist = hmeDist;
 638
 639 //                     Pout<< "Merging slave point; SECOND TRY "
 640 //                         << slavePatch.meshPoints()[pointI] << " local "
 641 //                         << pointI << " at "
 642 //                         << slaveLocalPoints[pointI] << " with master "
 643 //                         << masterPatch.meshPoints()[mergePoint] << " at "
 644 //                         << masterLocalPoints[mergePoint]
 645 //                         << ". hmeDist: " << mergeDist
 646 //                         << " mergeTol: " << mergeTol << endl;
 647                     }
 648                 }
 649
 650                 if (mergePoint > -1)
 651                 {
 652                     // Point is to be merged with master point
 653                     slavePointPointHits[pointI] = mergePoint;
 654                     curPoint = masterLocalPoints[mergePoint];
 655                     masterPointPointHits[mergePoint] = pointI;
 656
 657                     slavePointFaceHits[pointI] =
 658                         objectHit(true, slavePointFaceHits[pointI].hitObject());
 659
 660
 661                     // Disable edge merge
 662                     slavePointEdgeHits[pointI] = -1;
 663                 }
 664             }
 665         }
 666     }
 667
 668     if (debug)
 669     {
 670         Pout<< "Number of merged master points: " << nMergedPoints << nl
 671             << "Number of adjusted slave points: " << nMovedPoints << endl;
 672
 673         // Check for zero-length edges in slave projection
 674         scalar minEdgeLength = GREAT;
 675         scalar el = 0;
 676
 677         forAll (slaveEdges, edgeI)
 678         {
 679             el = slaveEdges[edgeI].mag(projectedSlavePoints);
 680
 681             if (el < SMALL)
 682             {
 683                 Pout<< "Point projection problems for edge: "
 684                     << slaveEdges[edgeI] << ". Length = " << el
 685                     << endl;
 686             }
 687
 688             minEdgeLength = min(minEdgeLength, el);
 689         }
 690
 691         if (minEdgeLength < SMALL)
 692         {
 693             FatalErrorIn
 694             (
 695                 "bool slidingInterface::projectPoints() const"
 696             )   << "Problem in point projection.  Short projected edge "
 697             << " after correction for object " << name()
 698             << abort(FatalError);
 699         }
 700     }
 701
 702 //     Pout << "slavePointEdgeHits: " << slavePointEdgeHits << endl;
 703
 704     // Insert the master points into closest slave edge if appropriate
 705
 706     // Algorithm:
 707     //    The face potentially interacts with all the points of the
 708     //    faces covering the path between its two ends.  Since we are
 709     //    examining an arbitrary shadow of the edge on a non-Euclidian
 710     //    surface, it is typically quite hard to do a geometric point
 711     //    search (under a shadow of a straight line).  Therefore, the
 712     //    search will be done topologically.
 713     //
 714     // I) Point collection
 715     // -------------------
 716     // 1) Grab the master faces to which the end points of the edge
 717     //    have been projected.
 718     // 2) Starting from the face containing the edge start, grab all
 719     //    its points and put them into the point lookup map.  Put the
 720     //    face onto the face lookup map.
 721     // 3) If the face of the end point is on the face lookup, complete
 722     //    the point collection step (this is checked during insertion.
 723     // 4) Start a new round of insertion.  Visit all faces in the face
 724     //    lookup and add their neighbours which are not already on the
 725     //    map.  Every time the new neighbour is found, add its points to
 726     //    the point lookup.  If the face of the end point is inserted,
 727     //    continue with the current roundof insertion and stop the
 728     //    algorithm.
 729     //
 730     // II) Point check
 731     // ---------------
 732     //    Grab all the points from the point collection and check them
 733     //    against the current edge.  If the edge-to-point distance is
 734     //    smaller than the smallest tolerance in the game (min of
 735     //    master point tolerance and left and right slave face around
 736     //    the edge tolerance) and the nearest point hits the edge, the
 737     //    master point will break the slave edge.  Check the actual
 738     //    distance of the original master position from the edge.  If
 739     //    the distance is smaller than a fraction of slave edge
 740     //    length, the hit is considered valid.  Store the slave edge
 741     //    index for every master point.
 742
 743     labelList masterPointEdgeHits(masterLocalPoints.size(), -1);
 744     scalarField masterPointEdgeDist(masterLocalPoints.size(), GREAT);
 745
 746     // Note.  "Processing slave edges" code is repeated twice in the
 747     // sliding intergace coupling in order to allow the point
 748     // projection to be done separately from the actual cutting.
 749     // Please change consistently with coupleSlidingInterface.C
 750     //
 751
 752     if (debug)
 753     {
 754         Pout << "Processing slave edges " << endl;
 755     }
 756
 757     // Create a map of faces the edge can interact with
 758     labelHashSet curFaceMap
 759     (
 760         nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
 761     );
 762
 763     labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
 764
 765     forAll (slaveEdges, edgeI)
 766     {
 767         const edge& curEdge = slaveEdges[edgeI];
 768
 769         //HJ: check for all edges even if both ends have missed
 770         //    Experimental.
 771 //         if
 772 //         (
 773 //             slavePointFaceHits[curEdge.start()].hit()
 774 //          || slavePointFaceHits[curEdge.end()].hit()
 775 //         )
 776         {
 777             // Clear the maps
 778             curFaceMap.clear();
 779             addedFaces.clear();
 780
 781             // Grab the faces for start and end points
 782             const label startFace =
 783                 slavePointFaceHits[curEdge.start()].hitObject();
 784             const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
 785
 786             // Insert the start face into the list
 787             curFaceMap.insert(startFace);
 788             addedFaces.insert(startFace);
 789
 790 //             Pout << "Doing edge " << edgeI << " or " << curEdge << " start: " << slavePointFaceHits[curEdge.start()].hitObject() << " end " << slavePointFaceHits[curEdge.end()].hitObject() << endl;
 791             // If the end face is on the list, the face collection is finished
 792             label nSweeps = 0;
 793             bool completed = false;
 794
 795             while (nSweeps < edgeFaceEscapeLimit_)
 796             {
 797                 nSweeps++;
 798
 799                 if (addedFaces.found(endFace))
 800                 {
 801                     completed = true;
 802                 }
 803
 804                 // Add all face neighbours of face in the map
 805                 const labelList cf = addedFaces.toc();
 806                 addedFaces.clear();
 807
 808                 forAll (cf, cfI)
 809                 {
 810                     const labelList& curNbrs = masterFaceFaces[cf[cfI]];
 811
 812                     forAll (curNbrs, nbrI)
 813                     {
 814                         if (!curFaceMap.found(curNbrs[nbrI]))
 815                         {
 816                             curFaceMap.insert(curNbrs[nbrI]);
 817                             addedFaces.insert(curNbrs[nbrI]);
 818                         }
 819                     }
 820                 }
 821
 822                 if (completed) break;
 823
 824                 if (debug)
 825                 {
 826                     Pout << ".";
 827                 }
 828             }
 829
 830             if (!completed)
 831             {
 832                 if (debug)
 833                 {
 834                     Pout << "x";
 835                 }
 836
 837                 // It is impossible to reach the end from the start, probably
 838                 // due to disconnected domain.  Do search in opposite direction
 839
 840                 label nReverseSweeps = 0;
 841
 842                 addedFaces.clear();
 843                 curFaceMap.insert(endFace);
 844                 addedFaces.insert(endFace);
 845
 846                 while (nReverseSweeps < edgeFaceEscapeLimit_)
 847                 {
 848                     nReverseSweeps++;
 849
 850                     if (addedFaces.found(startFace))
 851                     {
 852                         completed = true;
 853                     }
 854
 855                     // Add all face neighbours of face in the map
 856                     const labelList cf = addedFaces.toc();
 857                     addedFaces.clear();
 858
 859                     forAll (cf, cfI)
 860                     {
 861                         const labelList& curNbrs = masterFaceFaces[cf[cfI]];
 862
 863                         forAll (curNbrs, nbrI)
 864                         {
 865                             if (!curFaceMap.found(curNbrs[nbrI]))
 866                             {
 867                                 curFaceMap.insert(curNbrs[nbrI]);
 868                                 addedFaces.insert(curNbrs[nbrI]);
 869                             }
 870                         }
 871                     }
 872
 873                     if (completed) break;
 874
 875                     if (debug)
 876                     {
 877                         Pout << ".";
 878                     }
 879                 }
 880             }
 881
 882             if (completed)
 883             {
 884                 if (debug)
 885                 {
 886                     Pout << "+ ";
 887                 }
 888             }
 889             else
 890             {
 891                 if (debug)
 892                 {
 893                     Pout << "z ";
 894                 }
 895             }
 896
 897             // Collect the points
 898
 899             // Create a map of points the edge can interact with
 900             labelHashSet curPointMap
 901             (
 902                 nFacesPerSlaveEdge_*primitiveMesh::pointsPerFace_
 903             );
 904
 905             const labelList curFaces = curFaceMap.toc();
 906 //             Pout << "curFaces: " << curFaces << endl;
 907             forAll (curFaces, faceI)
 908             {
 909                 const face& f = masterLocalFaces[curFaces[faceI]];
 910
 911                 forAll (f, pointI)
 912                 {
 913                     curPointMap.insert(f[pointI]);
 914                 }
 915             }
 916
 917             const labelList curMasterPoints = curPointMap.toc();
 918
 919             // Check all the points against the edge.
 920
 921             linePointRef edgeLine = curEdge.line(projectedSlavePoints);
 922
 923             const vector edgeVec = edgeLine.vec();
 924             const scalar edgeMag = edgeLine.mag();
 925
 926             // Calculate actual distance involved in projection.  This
 927             // is used to reject master points out of reach.
 928             // Calculated as a combination of travel distance in projection and
 929             // edge length
 930             scalar slaveCatchDist =
 931                 edgeMasterCatchFraction_*edgeMag
 932               + 0.5*
 933                 (
 934                     mag
 935                     (
 936                         projectedSlavePoints[curEdge.start()]
 937                       - slaveLocalPoints[curEdge.start()]
 938                     )
 939                   + mag
 940                     (
 941                         projectedSlavePoints[curEdge.end()]
 942                       - slaveLocalPoints[curEdge.end()]
 943                     )
 944                 );
 945
 946             // The point merge distance needs to be measured in the
 947             // plane of the slave edge.  The unit vector is calculated
 948             // as a cross product of the edge vector and the edge
 949             // projection direction.  When checking for the distance
 950             // in plane, a minimum of the master-to-edge and
 951             // projected-master-to-edge distance is used, to avoid
 952             // problems with badly defined master planes.  HJ,
 953             // 17/Oct/2004
 954             vector edgeNormalInPlane =
 955                 edgeVec
 956               ^ (
 957                     slavePointNormals[curEdge.start()]
 958                   + slavePointNormals[curEdge.end()]
 959                 );
 960
 961             edgeNormalInPlane /= mag(edgeNormalInPlane);
 962
 963             forAll (curMasterPoints, pointI)
 964             {
 965                 const label cmp = curMasterPoints[pointI];
 966
 967                 // Skip the current point if the edge start or end has
 968                 // been adjusted onto in
 969                 if
 970                 (
 971                     slavePointPointHits[curEdge.start()] == cmp
 972                  || slavePointPointHits[curEdge.end()] == cmp
 973                  || masterPointPointHits[cmp] > -1
 974                 )
 975                 {
 976 // Pout << "Edge already snapped to point.  Skipping." << endl;
 977                     continue;
 978                 }
 979
 980                 // Check if the point actually hits the edge within bounds
 981                 pointHit edgeLineHit =
 982                     edgeLine.nearestDist(masterLocalPoints[cmp]);
 983
 984                 if (edgeLineHit.hit())
 985                 {
 986                     // If the distance to the line is smaller than
 987                     // the tolerance the master point needs to be
 988                     // inserted into the edge
 989
 990                     // Strict checking of slave cut to avoid capturing
 991                     // end points.
 992                     scalar cutOnSlave =
 993                         ((edgeLineHit.hitPoint() - edgeLine.start()) & edgeVec)
 994                         /sqr(edgeMag);
 995
 996                     scalar distInEdgePlane =
 997                         min
 998                         (
 999                             edgeLineHit.distance(),
1000                             mag
1001                             (
1002                                 (
1003                                     masterLocalPoints[cmp]
1004                                   - edgeLineHit.hitPoint()
1005                                 )
1006                               & edgeNormalInPlane
1007                             )
1008                         );
1009 //                     Pout << "master point: " << cmp
1010 //                         << " cutOnSlave " << cutOnSlave
1011 //                         << " distInEdgePlane: " << distInEdgePlane
1012 //                         << " tol1: " << pointMergeTol_*edgeMag
1013 //                         << " hitDist: " << edgeLineHit.distance()
1014 //                         << " tol2: " <<
1015 //                         min
1016 //                         (
1017 //                             slaveCatchDist,
1018 //                             masterPointEdgeDist[cmp]
1019 //                         ) << endl;
1020
1021                     // Not a point hit, check for edge
1022                     if
1023                     (
1024                         cutOnSlave > edgeEndCutoffTol_
1025                      && cutOnSlave < 1.0 - edgeEndCutoffTol_ // check edge cut
1026                      && distInEdgePlane < edgeMergeTol_*edgeMag // merge plane
1027                      && edgeLineHit.distance()
1028                       < min
1029                         (
1030                             slaveCatchDist,
1031                             masterPointEdgeDist[cmp]
1032                         )
1033                     )
1034                     {
1035                         if (debug)
1036                         {
1037                             if (masterPointEdgeHits[cmp] == -1)
1038                             {
1039                                 // First hit
1040                                 Pout << "m";
1041                             }
1042                             else
1043                             {
1044                                 // Repeat hit
1045                                 Pout << "M";
1046                             }
1047                         }
1048
1049                         // Snap to point onto edge
1050                         masterPointEdgeHits[cmp] = edgeI;
1051                         masterPointEdgeDist[cmp] = edgeLineHit.distance();
1052
1053 //                         Pout<< "Inserting master point "
1054 //                             << masterPatch.meshPoints()[cmp]
1055 //                             << " (" << cmp
1056 //                             << ") at " << masterLocalPoints[cmp]
1057 //                             << " into slave edge " << edgeI
1058 //                             << " " << curEdge
1059 //                             << " cutOnSlave: " << cutOnSlave
1060 //                             << " distInEdgePlane: " << distInEdgePlane
1061 //                             << ". dist: " << edgeLineHit.distance()
1062 //                             << " mergeTol: "
1063 //                             << edgeMergeTol_*edgeMag
1064 //                             << " other tol: " <<
1065 //                             min
1066 //                             (
1067 //                                 slaveCatchDist,
1068 //                                 masterPointEdgeDist[cmp]
1069 //                             )
1070 //                             << endl;
1071                     }
1072                 }
1073             }
1074         } // End if both ends missing
1075     } // End all slave edges
1076
1077     if (debug)
1078     {
1079         Pout << endl;
1080     }
1081 //     Pout << "masterPointEdgeHits: " << masterPointEdgeHits << endl;
1082
1083     if (debug)
1084     {
1085         Pout<< "bool slidingInterface::projectPoints() for object "
1086             << name() << " : "
1087             << "Finished projecting  points. Topology = ";
1088     }
1089
1090 //     Pout<< "slavePointPointHits: " << slavePointPointHits << nl
1091 //         << "slavePointEdgeHits: " << slavePointEdgeHits << nl
1092 //         << "slavePointFaceHits: " << slavePointFaceHits << nl
1093 //         << "masterPointEdgeHits: " << masterPointEdgeHits << endl;
1094
1095     // The four lists:
1096     // - slavePointPointHits
1097     // - slavePointEdgeHits
1098     // - slavePointFaceHits
1099     // - masterPointEdgeHits
1100     //   define how the two patches will be merged topologically.
1101     //   If the lists have not changed since the last merge, the
1102     //   sliding interface changes only geometrically and simple mesh
1103     //   motion will suffice.  Otherwise, a topological change is
1104     //   required.
1105
1106     // Compare the result with the current state
1107     if (!attached_)
1108     {
1109         // The mesh needs to change topologically
1110         trigger_ = true;
1111
1112         // Store the addressing arrays and projected points
1113         deleteDemandDrivenData(slavePointPointHitsPtr_);
1114         slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1115
1116         deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1117         slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1118
1119         deleteDemandDrivenData(slavePointFaceHitsPtr_);
1120         slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1121
1122         deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1123         masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1124
1125         if (debug)
1126         {
1127             Pout << "(Detached interface) changing." << endl;
1128         }
1129     }
1130     else
1131     {
1132         // Compare the lists against the stored lists.  If any of them
1133         // has changed, topological change will be executed.
1134         trigger_ = false;
1135
1136         if
1137         (
1138             !slavePointPointHitsPtr_
1139          || !slavePointEdgeHitsPtr_
1140          || !slavePointFaceHitsPtr_
1141          || !masterPointEdgeHitsPtr_
1142         )
1143         {
1144             // Must be restart.  Force topological change.
1145             deleteDemandDrivenData(slavePointPointHitsPtr_);
1146             slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1147
1148             deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1149             slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1150
1151             deleteDemandDrivenData(slavePointFaceHitsPtr_);
1152             slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1153
1154             deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1155             masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1156
1157             if (debug)
1158             {
1159                 Pout << "(Attached interface restart) changing." << endl;
1160             }
1161
1162             trigger_ = true;
1163             return trigger_;
1164         }
1165
1166         if (slavePointPointHits != (*slavePointPointHitsPtr_))
1167         {
1168             if (debug)
1169             {
1170                 Pout << "(Point projection) ";
1171             }
1172
1173             trigger_ = true;
1174         }
1175
1176         if (slavePointEdgeHits != (*slavePointEdgeHitsPtr_))
1177         {
1178             if (debug)
1179             {
1180                 Pout << "(Edge projection) ";
1181             }
1182
1183             trigger_ = true;
1184         }
1185
1186         // Comparison for face hits needs to exclude points that have hit
1187         // another point or edge
1188         bool faceHitsDifferent = false;
1189
1190         const List<objectHit>& oldPointFaceHits = *slavePointFaceHitsPtr_;
1191
1192         forAll (slavePointFaceHits, pointI)
1193         {
1194             if
1195             (
1196                 slavePointPointHits[pointI] < 0
1197              && slavePointEdgeHits[pointI] < 0
1198             )
1199             {
1200                 // This is a straight face hit
1201                 if (slavePointFaceHits[pointI] != oldPointFaceHits[pointI])
1202                 {
1203                     // Two lists are different
1204                     faceHitsDifferent = true;
1205                     break;
1206                 }
1207             }
1208         }
1209
1210         if (faceHitsDifferent)
1211         {
1212             if (debug)
1213             {
1214                 Pout << "(Face projection) ";
1215             }
1216
1217             trigger_ = true;
1218
1219         }
1220
1221         if (masterPointEdgeHits != (*masterPointEdgeHitsPtr_))
1222         {
1223             if (debug)
1224             {
1225                 Pout << "(Master point projection) ";
1226             }
1227
1228             trigger_ = true;
1229         }
1230
1231
1232         if (trigger_)
1233         {
1234             // Grab new data
1235             deleteDemandDrivenData(slavePointPointHitsPtr_);
1236             slavePointPointHitsPtr_ = new labelList(slavePointPointHits);
1237
1238             deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1239             slavePointEdgeHitsPtr_ = new labelList(slavePointEdgeHits);
1240
1241             deleteDemandDrivenData(slavePointFaceHitsPtr_);
1242             slavePointFaceHitsPtr_ = new List<objectHit>(slavePointFaceHits);
1243
1244             deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1245             masterPointEdgeHitsPtr_ = new labelList(masterPointEdgeHits);
1246
1247             if (debug)
1248             {
1249                 Pout << "changing." << endl;
1250             }
1251         }
1252         else
1253         {
1254             if (debug)
1255             {
1256                 Pout << "preserved." << endl;
1257             }
1258         }
1259     }
1260
1261     return trigger_;
1262 }
1263
1264
1265 void Foam::slidingInterface::clearPointProjection() const
1266 {
1267     deleteDemandDrivenData(slavePointPointHitsPtr_);
1268     deleteDemandDrivenData(slavePointEdgeHitsPtr_);
1269     deleteDemandDrivenData(slavePointFaceHitsPtr_);
1270     deleteDemandDrivenData(masterPointEdgeHitsPtr_);
1271
1272     deleteDemandDrivenData(projectedSlavePointsPtr_);
1273 }
1274
1275
1276 // ************************************************************************* //