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