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