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