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