src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 \*---------------------------------------------------------------------------*/
  25
  26 #include "face.H"
  27 #include "pointHit.H"
  28 #include "triPointRef.H"
  29 #include "line.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 // Return potential intersection with face with a ray starting
  34 // at p, direction n (does not need to be normalized)
  35 // Does face-center decomposition and returns triangle intersection
  36 // point closest to p.
  37
  38 // In case of miss the point is the nearest point intersection of the
  39 // face plane and the ray and the distance is the distance between the
  40 // intersection point and the nearest point on the face
  41
  42 Foam::pointHit Foam::face::ray
  43 (
  44     const point& p,
  45     const vector& n,
  46     const pointField& meshPoints,
  47     const intersection::algorithm alg,
  48     const intersection::direction dir
  49 ) const
  50 {
  51     // If the face is a triangle, do a direct calculation
  52     if (size() == 3)
  53     {
  54         return triPointRef
  55         (
  56             meshPoints[operator[](0)],
  57             meshPoints[operator[](1)],
  58             meshPoints[operator[](2)]
  59         ).ray(p, n, alg, dir);
  60     }
  61
  62     point ctr = Foam::average(points(meshPoints));
  63
  64     scalar nearestHitDist = GREAT;
  65
  66     scalar nearestMissDist = GREAT;
  67     bool eligible = false;
  68
  69     // Initialize to miss, distance = GREAT
  70     pointHit nearest(p);
  71
  72     const labelList& f = *this;
  73
  74     label nPoints = size();
  75
  76     point nextPoint = ctr;
  77
  78     for (label pI = 0; pI < nPoints; pI++)
  79     {
  80         nextPoint = meshPoints[f[fcIndex(pI)]];
  81
  82         // Note: for best accuracy, centre point always comes last
  83         //
  84         pointHit curHit = triPointRef
  85         (
  86             meshPoints[f[pI]],
  87             nextPoint,
  88             ctr
  89         ).ray(p, n, alg, dir);
  90
  91         if (curHit.hit())
  92         {
  93             if (Foam::mag(curHit.distance()) < Foam::mag(nearestHitDist))
  94             {
  95                 nearestHitDist = curHit.distance();
  96                 nearest.setHit();
  97                 nearest.setPoint(curHit.hitPoint());
  98             }
  99         }
 100         else if (!nearest.hit())
 101         {
 102             // Miss and no hit yet. Update miss statistics.
 103             if (curHit.eligibleMiss())
 104             {
 105                 eligible = true;
 106
 107                 // Miss distance is the distance between the plane intersection
 108                 // point and the nearest point of the triangle
 109                 scalar missDist =
 110                     Foam::mag
 111                     (
 112                         p + curHit.distance()*n
 113                       - curHit.missPoint()
 114                     );
 115
 116                 if (missDist < nearestMissDist)
 117                 {
 118                     nearestMissDist = missDist;
 119                     nearest.setDistance(curHit.distance());
 120                     nearest.setPoint(curHit.missPoint());
 121                 }
 122             }
 123         }
 124     }
 125
 126     if (nearest.hit())
 127     {
 128         nearest.setDistance(nearestHitDist);
 129     }
 130     else
 131     {
 132         // Haven't hit a single face triangle
 133         nearest.setMiss(eligible);
 134     }
 135
 136     return nearest;
 137 }
 138
 139
 140 Foam::pointHit Foam::face::intersection
 141 (
 142     const point& p,
 143     const vector& q,
 144     const point& ctr,
 145     const pointField& meshPoints,
 146     const intersection::algorithm alg,
 147     const scalar tol
 148 ) const
 149 {
 150     // If the face is a triangle, do a direct calculation
 151     if (size() == 3)
 152     {
 153         return triPointRef
 154         (
 155             meshPoints[operator[](0)],
 156             meshPoints[operator[](1)],
 157             meshPoints[operator[](2)]
 158         ).intersection(p, q, alg, tol);
 159     }
 160
 161     scalar nearestHitDist = VGREAT;
 162
 163     // Initialize to miss, distance = GREAT
 164     pointHit nearest(p);
 165
 166     const labelList& f = *this;
 167
 168     forAll(f, pI)
 169     {
 170         // Note: for best accuracy, centre point always comes last
 171         pointHit curHit = triPointRef
 172         (
 173             meshPoints[f[pI]],
 174             meshPoints[f[fcIndex(pI)]],
 175             ctr
 176         ).intersection(p, q, alg, tol);
 177
 178         if (curHit.hit())
 179         {
 180             if (Foam::mag(curHit.distance()) < Foam::mag(nearestHitDist))
 181             {
 182                 nearestHitDist = curHit.distance();
 183                 nearest.setHit();
 184                 nearest.setPoint(curHit.hitPoint());
 185             }
 186         }
 187     }
 188
 189     if (nearest.hit())
 190     {
 191         nearest.setDistance(nearestHitDist);
 192     }
 193
 194     return nearest;
 195 }
 196
 197
 198 Foam::pointHit Foam::face::nearestPoint
 199 (
 200     const point& p,
 201     const pointField& meshPoints
 202 ) const
 203 {
 204     // Dummy labels
 205     label nearType = -1;
 206     label nearLabel = -1;
 207
 208     return nearestPointClassify(p, meshPoints, nearType, nearLabel);
 209 }
 210
 211
 212 Foam::pointHit Foam::face::nearestPointClassify
 213 (
 214     const point& p,
 215     const pointField& meshPoints,
 216     label& nearType,
 217     label& nearLabel
 218 ) const
 219 {
 220     // If the face is a triangle, do a direct calculation
 221     if (size() == 3)
 222     {
 223         return triPointRef
 224         (
 225             meshPoints[operator[](0)],
 226             meshPoints[operator[](1)],
 227             meshPoints[operator[](2)]
 228         ).nearestPointClassify(p, nearType, nearLabel);
 229     }
 230
 231     const face& f = *this;
 232     point ctr = centre(meshPoints);
 233
 234     // Initialize to miss, distance=GREAT
 235     pointHit nearest(p);
 236
 237     nearType = -1;
 238     nearLabel = -1;
 239
 240     label nPoints = f.size();
 241
 242     point nextPoint = ctr;
 243
 244     for (label pI = 0; pI < nPoints; pI++)
 245     {
 246         nextPoint = meshPoints[f[fcIndex(pI)]];
 247
 248         label tmpNearType = -1;
 249         label tmpNearLabel = -1;
 250
 251         // Note: for best accuracy, centre point always comes last
 252         triPointRef tri
 253         (
 254             meshPoints[f[pI]],
 255             nextPoint,
 256             ctr
 257         );
 258
 259         pointHit curHit = tri.nearestPointClassify
 260         (
 261             p,
 262             tmpNearType,
 263             tmpNearLabel
 264         );
 265
 266         if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
 267         {
 268             nearest.setDistance(curHit.distance());
 269
 270             // Assume at first that the near type is NONE on the
 271             // triangle (i.e. on the face of the triangle) then it is
 272             // therefore also for the face.
 273
 274             nearType = NONE;
 275
 276             if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
 277             {
 278                 // If the triangle edge label is 0, then this is also
 279                 // an edge of the face, if not, it is on the face
 280
 281                 nearType = EDGE;
 282
 283                 nearLabel = pI;
 284             }
 285             else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
 286             {
 287                 // If the triangle point label is 0 or 1, then this is
 288                 // also a point of the face, if not, it is on the face
 289
 290                 nearType = POINT;
 291
 292                 nearLabel = pI + tmpNearLabel;
 293             }
 294
 295             if (curHit.hit())
 296             {
 297                 nearest.setHit();
 298                 nearest.setPoint(curHit.hitPoint());
 299             }
 300             else
 301             {
 302                 // In nearest point, miss is always eligible
 303                 nearest.setMiss(true);
 304                 nearest.setPoint(curHit.missPoint());
 305             }
 306         }
 307     }
 308
 309     return nearest;
 310 }
 311
 312
 313 // ************************************************************************* //