src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C

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