src/dynamicMesh/polyMeshModifiers/slidingInterface/slidingInterfaceProjectPoints.C

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