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