src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchProjectPoints.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     For every point on the patch find the closest face on the target side.
  26     Return a target face label for each patch point
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "boolList.H"
  31 #include "PointHitTemplate.H"
  32 #include "objectHit.H"
  33 #include "bandCompression.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 template
  38 <
  39     class Face,
  40     template<class> class FaceList,
  41     class PointField,
  42     class PointType
  43 >
  44 template <class ToPatch>
  45 Foam::List<Foam::objectHit>
  46 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
  47 projectPoints
  48 (
  49     const ToPatch& targetPatch,
  50     const Field<PointType>& projectionDirection,
  51     const intersection::algorithm alg,
  52     const intersection::direction dir
  53 ) const
  54 {
  55     // The current patch is slave, i.e. it is being projected onto the target
  56
  57     if (projectionDirection.size() != nPoints())
  58     {
  59         FatalErrorIn
  60         (
  61             "PrimitivePatch<Face, FaceList, PointField, PointType>::"
  62             "projectPoints(const PrimitivePatch& "
  63             ", const Field<PointType>&) const"
  64         )   << "Projection direction field does not correspond to "
  65             << "patch points." << endl
  66             << "Size: " << projectionDirection.size()
  67             << " Number of points: " << nPoints()
  68             << abort(FatalError);
  69     }
  70
  71     const labelList& slavePointOrder = localPointOrder();
  72
  73     const labelList& slaveMeshPoints = meshPoints();
  74
  75     // Result
  76     List<objectHit> result(nPoints());
  77
  78     // If there are no faces in the master patch, return: all misses
  79     if (targetPatch.size() == 0)
  80     {
  81         // Null-constructed object hit is a miss
  82         return result;
  83     }
  84
  85     const labelListList& masterFaceFaces = targetPatch.faceFaces();
  86
  87     const ToPatch& masterFaces = targetPatch;
  88
  89     const Field<PointType>& masterPoints = targetPatch.points();
  90
  91     // Estimate face centre of target side
  92     Field<PointType> masterFaceCentres(targetPatch.size());
  93
  94     forAll (masterFaceCentres, faceI)
  95     {
  96         masterFaceCentres[faceI] =
  97             average(masterFaces[faceI].points(masterPoints));
  98     }
  99
 100     // Algorithm:
 101     // Loop through all points of the slave side. For every point find the
 102     // radius for the current contact face. If the contact point falls inside
 103     // the face and the radius is smaller than for all neighbouring faces,
 104     // the contact is found. If not, visit the neighbour closest to the
 105     // calculated contact point. If a single master face is visited more than
 106     // twice, initiate n-squared search.
 107
 108     if (debug)
 109     {
 110         Info<< "N-squared projection set for all points" << endl;
 111     }
 112
 113     label curFace = 0;
 114     label nNSquaredSearches = 0;
 115
 116     forAll (slavePointOrder, pointI)
 117     {
 118         // Pick up slave point and direction
 119         const label curLocalPointLabel = slavePointOrder[pointI];
 120
 121         const PointType& curPoint =
 122             points_[slaveMeshPoints[curLocalPointLabel]];
 123
 124         const PointType& curProjectionDir =
 125             projectionDirection[curLocalPointLabel];
 126
 127         bool closer;
 128
 129         boolList visitedTargetFace(targetPatch.size(), false);
 130         bool doNSquaredSearch = false;
 131
 132         bool foundEligible = false;
 133
 134         scalar sqrDistance = GREAT;
 135
 136         // Force the full search for the first point to ensure good
 137         // starting face
 138         if (pointI == 0)
 139         {
 140             doNSquaredSearch = true;
 141         }
 142         else
 143         {
 144             do
 145             {
 146                 closer = false;
 147                 doNSquaredSearch = false;
 148
 149                 // Calculate intersection with curFace
 150                 PointHit<PointType> curHit =
 151                     masterFaces[curFace].ray
 152                     (
 153                         curPoint,
 154                         curProjectionDir,
 155                         masterPoints,
 156                         alg,
 157                         dir
 158                     );
 159
 160                 visitedTargetFace[curFace] = true;
 161
 162                 if (curHit.hit())
 163                 {
 164                     result[curLocalPointLabel] = objectHit(true, curFace);
 165
 166                     break;
 167                 }
 168                 else
 169                 {
 170                     // If a new miss is eligible, it is closer than
 171                     // any previous eligible miss (due to surface walk)
 172
 173                     // Only grab the miss if it is eligible
 174                     if (curHit.eligibleMiss())
 175                     {
 176                         foundEligible = true;
 177                         result[curLocalPointLabel] = objectHit(false, curFace);
 178                     }
 179
 180                     // Find the next likely face for intersection
 181
 182                     // Calculate the miss point on the plane of the
 183                     // face.  This is cooked (illogical!) for fastest
 184                     // surface walk.
 185                     //
 186                     PointType missPlanePoint =
 187                         curPoint + curProjectionDir*curHit.distance();
 188
 189                     const labelList& masterNbrs = masterFaceFaces[curFace];
 190
 191                     sqrDistance =
 192                         magSqr(missPlanePoint - masterFaceCentres[curFace]);
 193
 194                     forAll (masterNbrs, nbrI)
 195                     {
 196                         if
 197                         (
 198                             magSqr
 199                             (
 200                                 missPlanePoint
 201                               - masterFaceCentres[masterNbrs[nbrI]]
 202                             )
 203                          <= sqrDistance
 204                         )
 205                         {
 206                             closer = true;
 207                             curFace = masterNbrs[nbrI];
 208                         }
 209                     }
 210
 211                     if (visitedTargetFace[curFace])
 212                     {
 213                         // This face has already been visited.
 214                         // Execute n-squared search
 215                         doNSquaredSearch = true;
 216                         break;
 217                     }
 218                 }
 219
 220                 if (debug) Info<< ".";
 221             } while (closer);
 222         }
 223
 224         if (doNSquaredSearch || !foundEligible)
 225         {
 226             nNSquaredSearches++;
 227
 228             if (debug)
 229             {
 230                 Info<< "p " << curLocalPointLabel << ": ";
 231             }
 232
 233             // Improved n-sqared search algorithm.  HJ, 25/Nov/2006
 234             result[curLocalPointLabel] = objectHit(false, -1);
 235             scalar minMissDistance = GREAT;
 236             scalar minHitDistance = GREAT;
 237             bool hitFound = false;
 238
 239             forAll (masterFaces, faceI)
 240             {
 241                 PointHit<PointType> curHit =
 242                     masterFaces[faceI].ray
 243                     (
 244                         curPoint,
 245                         curProjectionDir,
 246                         masterPoints,
 247                         alg,
 248                         dir
 249                     );
 250
 251                 if (curHit.hit())
 252                 {
 253                     // Calculate min distance
 254                     scalar hitDist =
 255                         Foam::mag(curHit.hitPoint() - curPoint);
 256
 257                     if (hitDist < minHitDistance)
 258                     {
 259                         result[curLocalPointLabel] = objectHit(true, faceI);
 260
 261                         hitFound = true;
 262                         curFace = faceI;
 263                         minHitDistance = hitDist;
 264                     }
 265                 }
 266                 else if (curHit.eligibleMiss() && !hitFound)
 267                 {
 268                     // Calculate min distance
 269                     scalar missDist =
 270                         Foam::mag(curHit.missPoint() - curPoint);
 271
 272                     if (missDist < minMissDistance)
 273                     {
 274                         minMissDistance = missDist;
 275
 276                         result[curLocalPointLabel] = objectHit(false, faceI);
 277                         curFace = faceI;
 278                     }
 279                 }
 280             }
 281
 282             if (debug)
 283             {
 284                 Info<< result[curLocalPointLabel]
 285                     << " hit = " << minHitDistance
 286                     << " miss = " << minMissDistance << nl;
 287             }
 288         }
 289         else
 290         {
 291             if (debug) Info<< "x";
 292         }
 293     }
 294
 295     if (debug)
 296     {
 297         Info<< nl << "Executed " << nNSquaredSearches
 298             << " n-squared searches out of total of "
 299             << nPoints() << endl;
 300     }
 301
 302     return result;
 303 }
 304
 305
 306 template
 307 <
 308     class Face,
 309     template<class> class FaceList,
 310     class PointField,
 311     class PointType
 312 >
 313 template <class ToPatch>
 314 Foam::List<Foam::objectHit>
 315 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 316 projectFaceCentres
 317 (
 318     const ToPatch& targetPatch,
 319     const Field<PointType>& projectionDirection,
 320     const intersection::algorithm alg,
 321     const intersection::direction dir
 322 ) const
 323 {
 324     // The current patch is slave, i.e. it is being projected onto the target
 325
 326     if (projectionDirection.size() != this->size())
 327     {
 328         FatalErrorIn
 329         (
 330             "labelList PrimitivePatch<Face, FaceList, PointField, PointType>::"
 331             "projectFaceCentres(const PrimitivePatch& "
 332             ", const Field<PointType>&) const"
 333         )   << "Projection direction field does not correspond to patch faces."
 334             << endl << "Size: " << projectionDirection.size()
 335             << " Number of points: " << this->size()
 336             << abort(FatalError);
 337     }
 338
 339     labelList slaveFaceOrder = bandCompression(faceFaces());
 340
 341     // calculate master face centres
 342     Field<PointType> masterFaceCentres(targetPatch.size());
 343
 344     const labelListList& masterFaceFaces = targetPatch.faceFaces();
 345
 346     const ToPatch& masterFaces = targetPatch;
 347
 348     const Field<PointType>& masterPoints = targetPatch.points();
 349
 350     forAll (masterFaceCentres, faceI)
 351     {
 352         masterFaceCentres[faceI] =
 353             masterFaces[faceI].centre(masterPoints);
 354     }
 355
 356     // Result
 357     List<objectHit> result(this->size());
 358
 359     // If there are no faces in the master patch, return: all misses
 360     if (targetPatch.size() == 0)
 361     {
 362         // Null-constructed object hit is a miss
 363         return result;
 364     }
 365
 366     const PrimitivePatch<Face, FaceList, PointField>& slaveFaces = *this;
 367     const vectorField& slaveGlobalPoints = points();
 368
 369     // Algorithm:
 370     // Loop through all points of the slave side. For every point find the
 371     // radius for the current contact face. If the contact point falls inside
 372     // the face and the radius is smaller than for all neighbouring faces,
 373     // the contact is found. If not, visit the neighbour closest to the
 374     // calculated contact point. If a single master face is visited more than
 375     // twice, initiate n-squared search.
 376
 377     label curFace = 0;
 378     label nNSquaredSearches = 0;
 379
 380     forAll (slaveFaceOrder, faceI)
 381     {
 382         // pick up slave point and direction
 383         const label curLocalFaceLabel = slaveFaceOrder[faceI];
 384
 385         const point& curFaceCentre =
 386             slaveFaces[curLocalFaceLabel].centre(slaveGlobalPoints);
 387
 388         const vector& curProjectionDir =
 389             projectionDirection[curLocalFaceLabel];
 390
 391         bool closer;
 392
 393         boolList visitedTargetFace(targetPatch.size(), false);
 394         bool doNSquaredSearch = false;
 395
 396         bool foundEligible = false;
 397
 398         scalar sqrDistance = GREAT;
 399
 400         // Force the full search for the first point to ensure good
 401         // starting face
 402         if (faceI == 0 || nSquaredProjection_() > 0)
 403         {
 404             doNSquaredSearch = true;
 405         }
 406         else
 407         {
 408             do
 409             {
 410                 closer = false;
 411                 doNSquaredSearch = false;
 412
 413                 // Calculate intersection with curFace
 414                 PointHit<PointType> curHit =
 415                     masterFaces[curFace].ray
 416                     (
 417                         curFaceCentre,
 418                         curProjectionDir,
 419                         masterPoints,
 420                         alg,
 421                         dir
 422                     );
 423
 424                 visitedTargetFace[curFace] = true;
 425
 426                 if (curHit.hit())
 427                 {
 428                     result[curLocalFaceLabel] = objectHit(true, curFace);
 429
 430                     break;
 431                 }
 432                 else
 433                 {
 434                     // If a new miss is eligible, it is closer than
 435                     // any previous eligible miss (due to surface walk)
 436
 437                     // Only grab the miss if it is eligible
 438                     if (curHit.eligibleMiss())
 439                     {
 440                         foundEligible = true;
 441                         result[curLocalFaceLabel] = objectHit(false, curFace);
 442                     }
 443
 444                     // Find the next likely face for intersection
 445
 446                     // Calculate the miss point.  This is
 447                     // cooked (illogical!) for fastest surface walk.
 448                     //
 449                     PointType missPlanePoint =
 450                         curFaceCentre + curProjectionDir*curHit.distance();
 451
 452                     sqrDistance =
 453                         magSqr(missPlanePoint - masterFaceCentres[curFace]);
 454
 455                     const labelList& masterNbrs = masterFaceFaces[curFace];
 456
 457                     forAll (masterNbrs, nbrI)
 458                     {
 459                         if
 460                         (
 461                             magSqr
 462                             (
 463                                 missPlanePoint
 464                               - masterFaceCentres[masterNbrs[nbrI]]
 465                             )
 466                          <= sqrDistance
 467                         )
 468                         {
 469                             closer = true;
 470                             curFace = masterNbrs[nbrI];
 471                         }
 472                     }
 473
 474                     if (visitedTargetFace[curFace])
 475                     {
 476                         // This face has already been visited.
 477                         // Execute n-squared search
 478                         doNSquaredSearch = true;
 479                         break;
 480                     }
 481                 }
 482
 483                 if (debug) Info<< ".";
 484             } while (closer);
 485         }
 486
 487         if (doNSquaredSearch || !foundEligible)
 488         {
 489             nNSquaredSearches++;
 490
 491             if (debug)
 492             {
 493                 Info<< "p " << curLocalFaceLabel << ": ";
 494             }
 495
 496             // Improved n-sqared search algorithm.  HJ, 25/Nov/2006
 497             result[curLocalFaceLabel] = objectHit(false, -1);
 498             scalar minMissDistance = GREAT;
 499             scalar minHitDistance = GREAT;
 500             bool hitFound = false;
 501
 502             forAll (masterFaces, faceI)
 503             {
 504                 PointHit<PointType> curHit =
 505                     masterFaces[faceI].ray
 506                     (
 507                         curFaceCentre,
 508                         curProjectionDir,
 509                         masterPoints,
 510                         alg,
 511                         dir
 512                     );
 513
 514                 if (curHit.hit())
 515                 {
 516                     // Calculate min distance
 517                     scalar hitDist =
 518                         Foam::mag(curHit.hitPoint() - curFaceCentre);
 519
 520                     if (hitDist < minHitDistance)
 521                     {
 522                         result[curLocalFaceLabel] = objectHit(true, faceI);
 523
 524                         hitFound = true;
 525                         curFace = faceI;
 526                         minHitDistance = hitDist;
 527                     }
 528                 }
 529                 else if (curHit.eligibleMiss() && !hitFound)
 530                 {
 531                     // Calculate min distance
 532                     scalar missDist =
 533                         Foam::mag(curHit.missPoint() - curFaceCentre);
 534
 535                     if (missDist < minMissDistance)
 536                     {
 537                         minMissDistance = missDist;
 538
 539                         result[curLocalFaceLabel] = objectHit(false, faceI);
 540                         curFace = faceI;
 541                     }
 542                 }
 543             }
 544
 545             if (debug)
 546             {
 547                 Info << "r = " << result[curLocalFaceLabel] << nl;
 548             }
 549         }
 550         else
 551         {
 552             if (debug) Info<< "x";
 553         }
 554     }
 555
 556     if (debug)
 557     {
 558         Info<< nl << "Executed " << nNSquaredSearches
 559             << " n-squared searches out of total of "
 560             << this->size() << endl;
 561     }
 562
 563     return result;
 564 }
 565
 566
 567 // ************************************************************************* //