src/OpenFOAM/meshes/primitiveShapes/line/line.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "line.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  35
  36 template<>
  37 scalar line<point2D, const point2D&>::nearestDist
  38 (
  39     const line<point2D, const point2D&>& e,
  40     point2D& thisPt,
  41     point2D& edgePt
  42 ) const
  43 {
  44     vector2D u = end()-start();
  45     vector2D v = e.end()-e.start();
  46     vector2D w = start()-e.start();
  47
  48     scalar d = u.perp(v);
  49
  50     if (Foam::mag(d) > VSMALL)
  51     {
  52         scalar s = v.perp(w) / d;
  53
  54         if (s <= SMALL)
  55         {
  56             thisPt = start();
  57         }
  58         else if (s >= (1-SMALL))
  59         {
  60             thisPt = end();
  61         }
  62         else
  63         {
  64             thisPt = start()+s*u;
  65         }
  66
  67
  68         scalar t = u.perp(w) / d;
  69
  70         if (t <= SMALL)
  71         {
  72             edgePt = e.start();
  73         }
  74         else if (t >= (1-SMALL))
  75         {
  76             edgePt = e.end();
  77         }
  78         else
  79         {
  80             edgePt = e.start()+t*v;
  81         }
  82     }
  83     else
  84     {
  85         // Parallel lines. Find overlap of both lines by projecting onto
  86         // direction vector (now equal for both lines).
  87
  88         scalar edge0 = e.start() & u;
  89         scalar edge1 = e.end() & u;
  90         bool edgeOrder = edge0 < edge1;
  91
  92         scalar minEdge = (edgeOrder ? edge0 : edge1);
  93         scalar maxEdge = (edgeOrder ? edge1 : edge0);
  94         const point2D& minEdgePt = (edgeOrder ? e.start() : e.end());
  95         const point2D& maxEdgePt = (edgeOrder ? e.end() : e.start());
  96
  97         scalar this0 = start() & u;
  98         scalar this1 = end() & u;
  99         bool thisOrder = this0 < this1;
 100
 101         scalar minThis = min(this0, this1);
 102         scalar maxThis = max(this1, this0);
 103         const point2D& minThisPt = (thisOrder ? start() : end());
 104         const point2D& maxThisPt = (thisOrder ? end() : start());
 105
 106         if (maxEdge < minThis)
 107         {
 108             // edge completely below *this
 109             edgePt = maxEdgePt;
 110             thisPt = minThisPt;
 111         }
 112         else if (maxEdge < maxThis)
 113         {
 114             // maxEdge inside interval of *this
 115             edgePt = maxEdgePt;
 116             thisPt = nearestDist(edgePt).rawPoint();
 117         }
 118         else
 119         {
 120             // maxEdge outside. Check if minEdge inside.
 121             if (minEdge < minThis)
 122             {
 123                 // Edge completely envelops this. Take any this point and
 124                 // determine nearest on edge.
 125                 thisPt = minThisPt;
 126                 edgePt = e.nearestDist(thisPt).rawPoint();
 127             }
 128             else if (minEdge < maxThis)
 129             {
 130                 // minEdge inside this interval.
 131                 edgePt = minEdgePt;
 132                 thisPt = nearestDist(edgePt).rawPoint();
 133             }
 134             else
 135             {
 136                 // minEdge outside this interval
 137                 edgePt = minEdgePt;
 138                 thisPt = maxThisPt;
 139             }
 140         }
 141     }
 142
 143     return Foam::mag(thisPt - edgePt);
 144 }
 145
 146
 147 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 148
 149 } // End namespace Foam
 150
 151 // ************************************************************************* //