src/dynamicMesh/dynamicTopoFvMesh/convexSetAlgorithm/triIntersectionI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Implemented by
  25     Sandeep Menon
  26     University of Massachusetts Amherst
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 namespace Foam
  31 {
  32
  33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  34
  35 // Compute clip-planes
  36 inline void triIntersection::computeClipPlanes()
  37 {
  38     // Face normal
  39     tNorm_ = (clipTri_[1] - clipTri_[0]) ^ (clipTri_[2] - clipTri_[0]);
  40
  41     // Normalize
  42     tNorm_ /= mag(tNorm_) + VSMALL;
  43
  44     // Edge normals
  45     clipPlanes_[0].first() = ((clipTri_[1] - clipTri_[0]) ^ tNorm_);
  46     clipPlanes_[1].first() = ((clipTri_[2] - clipTri_[1]) ^ tNorm_);
  47     clipPlanes_[2].first() = ((clipTri_[0] - clipTri_[2]) ^ tNorm_);
  48
  49     // Normalize
  50     clipPlanes_[0].first() /= mag(clipPlanes_[0].first()) + VSMALL;
  51     clipPlanes_[1].first() /= mag(clipPlanes_[1].first()) + VSMALL;
  52     clipPlanes_[2].first() /= mag(clipPlanes_[2].first()) + VSMALL;
  53
  54     // Determine plane constants
  55     clipPlanes_[0].second() = (clipTri_[0] & clipPlanes_[0].first());
  56     clipPlanes_[1].second() = (clipTri_[1] & clipPlanes_[1].first());
  57     clipPlanes_[2].second() = (clipTri_[2] & clipPlanes_[2].first());
  58 }
  59
  60
  61 // Split and decompose
  62 inline void triIntersection::splitAndDecompose
  63 (
  64     const label triPlaneIndex,
  65     FixedList<point, 3>& tri,
  66     DynamicList<FixedList<point, 3> >& decompTris
  67 ) const
  68 {
  69     FixedList<scalar, 3> C;
  70     FixedList<vector, 3> tmpTri;
  71     FixedList<label, 3> pos, neg, zero;
  72     label i = 0, nPos = 0, nNeg = 0, nZero = 0;
  73
  74     // Fetch reference to plane
  75     const hPlane& triPlane = clipPlanes_[triPlaneIndex];
  76
  77     for (i = 0; i < 3; ++i)
  78     {
  79         // Compute distance to plane
  80         C[i] = (tri[i] & triPlane.first()) - triPlane.second();
  81
  82         if (C[i] > 0.0)
  83         {
  84             pos[nPos++] = i;
  85         }
  86         else
  87         if (C[i] < 0.0)
  88         {
  89             neg[nNeg++] = i;
  90         }
  91         else
  92         {
  93             zero[nZero++] = i;
  94         }
  95     }
  96
  97     if (nNeg == 0)
  98     {
  99         return;
 100     }
 101
 102     if (nPos == 0)
 103     {
 104         decompTris.append(tri);
 105         return;
 106     }
 107
 108     // Triangle is split by plane.  Determine how it is split and how to
 109     // decompose the negative-side portion into triangles (3 cases).
 110     scalar w0, w1, invCDiff;
 111     vector intp[3];
 112
 113     if (nPos == 2)
 114     {
 115         // ++-
 116         for (i = 0; i < nPos; ++i)
 117         {
 118             invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
 119
 120             w0 = -C[neg[0]] * invCDiff;
 121             w1 = +C[pos[i]] * invCDiff;
 122
 123             tri[pos[i]] = (w0 * tri[pos[i]]) + (w1 * tri[neg[0]]);
 124         }
 125
 126         decompTris.append(tri);
 127     }
 128     else
 129     if (nPos == 1)
 130     {
 131         if (nNeg == 2)
 132         {
 133             // +--
 134             for (i = 0; i < nNeg; ++i)
 135             {
 136                 invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
 137
 138                 w0 = -C[neg[i]] * invCDiff;
 139                 w1 = +C[pos[0]] * invCDiff;
 140
 141                 intp[i] = (w0 * tri[pos[0]]) + (w1 * tri[neg[i]]);
 142             }
 143
 144             tri[pos[0]] = intp[0];
 145
 146             decompTris.append(tri);
 147
 148             tmpTri[0] = intp[0];
 149             tmpTri[1] = tri[neg[1]];
 150             tmpTri[2] = intp[1];
 151
 152             decompTris.append(tmpTri);
 153         }
 154         else
 155         {
 156             // +-0
 157             invCDiff = (1.0 / (C[pos[0]] - C[neg[0]]));
 158
 159             w0 = -C[neg[0]] * invCDiff;
 160             w1 = +C[pos[0]] * invCDiff;
 161
 162             tri[pos[0]] = (w0 * tri[pos[0]]) + (w1 * tri[neg[0]]);
 163
 164             decompTris.append(tri);
 165         }
 166     }
 167 }
 168
 169
 170 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 171
 172 inline triIntersection::triIntersection(const FixedList<point, 3>& clipTri)
 173 :
 174     clipTri_(clipTri),
 175     inside_(10),
 176     allTris_(10)
 177 {
 178     // Pre-compute clipping planes
 179     computeClipPlanes();
 180 }
 181
 182
 183 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 184
 185 inline triIntersection::~triIntersection()
 186 {}
 187
 188
 189 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 190
 191 // Evaluate for intersections
 192 inline bool triIntersection::evaluate(const FixedList<point, 3>& subjectTri)
 193 {
 194     // Clear lists
 195     inside_.clear();
 196     allTris_.clear();
 197
 198     // Project subject triangle to clipping triangle plane
 199     FixedList<point, 3> sTP(subjectTri);
 200
 201     for (label i = 0; i < 3; i++)
 202     {
 203         // Relative to zeroth face-point
 204         sTP[i] -= clipTri_[0];
 205
 206         // Project to face
 207         sTP[i] = clipTri_[0] + (sTP[i] - (sTP[i] & tNorm_)*tNorm_);
 208     }
 209
 210     // Add initial triangle to list
 211     allTris_.append(sTP);
 212
 213     // Clip second triangle against planes of clipping triangle
 214     for (label i = 0; i < 3; i++)
 215     {
 216         forAll(allTris_, triI)
 217         {
 218             splitAndDecompose(i, allTris_[triI], inside_);
 219         }
 220
 221         // Prep for next clipping plane
 222         allTris_ = inside_;
 223         inside_.clear();
 224     }
 225
 226     return (allTris_.size() > 0);
 227 }
 228
 229
 230 // Return intersections
 231 inline const DynamicList<FixedList<point, 3> >&
 232 triIntersection::getIntersection() const
 233 {
 234     return allTris_;
 235 }
 236
 237
 238 //- Evaluate and return area / centroid
 239 inline void triIntersection::getAreaAndCentre
 240 (
 241     scalar& area,
 242     vector& centre
 243 ) const
 244 {
 245     area = 0.0;
 246     centre = vector::zero;
 247
 248     forAll(allTris_, triI)
 249     {
 250         const FixedList<point, 3>& t = allTris_[triI];
 251
 252         // Calculate area (no check for orientation)
 253         scalar tA = Foam::mag(0.5 * ((t[1] - t[0]) ^ (t[2] - t[0])));
 254
 255         // Calculate centroid
 256         vector tC = (1.0 / 3.0) * (t[0] + t[1] + t[2]);
 257
 258         area += tA;
 259         centre += (tA * tC);
 260     }
 261
 262     centre /= area + VSMALL;
 263 }
 264
 265
 266 }
 267
 268 // ************************************************************************* //