src/dynamicMesh/dynamicTopoFvMesh/convexSetAlgorithm/tetIntersectionI.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 tetIntersection::computeClipPlanes()
  37 {
  38     // Define edge vectors
  39     vector edge10 = clipTet_[1] - clipTet_[0];
  40     vector edge20 = clipTet_[2] - clipTet_[0];
  41     vector edge30 = clipTet_[3] - clipTet_[0];
  42     vector edge21 = clipTet_[2] - clipTet_[1];
  43     vector edge31 = clipTet_[3] - clipTet_[1];
  44
  45     // Cross-products
  46     clipPlanes_[0].first() = (edge20 ^ edge10);
  47     clipPlanes_[1].first() = (edge10 ^ edge30);
  48     clipPlanes_[2].first() = (edge30 ^ edge20);
  49     clipPlanes_[3].first() = (edge21 ^ edge31);
  50
  51     // Normalize
  52     clipPlanes_[0].first() /= mag(clipPlanes_[0].first()) + VSMALL;
  53     clipPlanes_[1].first() /= mag(clipPlanes_[1].first()) + VSMALL;
  54     clipPlanes_[2].first() /= mag(clipPlanes_[2].first()) + VSMALL;
  55     clipPlanes_[3].first() /= mag(clipPlanes_[3].first()) + VSMALL;
  56
  57     // Compute magnitude of clipping tetrahedron
  58     clipTetMag_ = (1.0 / 6.0) * (edge10 & clipPlanes_[3].first());
  59
  60     if (clipTetMag_ < 0.0)
  61     {
  62         // Reverse normal directions
  63         clipPlanes_[0].first() = -clipPlanes_[0].first();
  64         clipPlanes_[1].first() = -clipPlanes_[1].first();
  65         clipPlanes_[2].first() = -clipPlanes_[2].first();
  66         clipPlanes_[3].first() = -clipPlanes_[3].first();
  67
  68         // Reverse sign
  69         clipTetMag_ = mag(clipTetMag_);
  70     }
  71
  72     // Determine plane constants
  73     clipPlanes_[0].second() = (clipTet_[0] & clipPlanes_[0].first());
  74     clipPlanes_[1].second() = (clipTet_[1] & clipPlanes_[1].first());
  75     clipPlanes_[2].second() = (clipTet_[2] & clipPlanes_[2].first());
  76     clipPlanes_[3].second() = (clipTet_[3] & clipPlanes_[3].first());
  77 }
  78
  79
  80 // Split and decompose
  81 inline void tetIntersection::splitAndDecompose
  82 (
  83     const label tetPlaneIndex,
  84     FixedList<point, 4>& tetra,
  85     DynamicList<FixedList<point, 4> >& decompTets
  86 ) const
  87 {
  88     FixedList<scalar, 4> C;
  89     FixedList<vector, 4> tmpTetra;
  90     FixedList<label, 4> pos, neg, zero;
  91     label i = 0, nPos = 0, nNeg = 0, nZero = 0;
  92
  93     // Fetch reference to plane
  94     const hPlane& tetPlane = clipPlanes_[tetPlaneIndex];
  95
  96     for (i = 0; i < 4; ++i)
  97     {
  98         // Compute distance to plane
  99         C[i] = (tetra[i] & tetPlane.first()) - tetPlane.second();
 100
 101         if (C[i] > 0.0)
 102         {
 103             pos[nPos++] = i;
 104         }
 105         else
 106         if (C[i] < 0.0)
 107         {
 108             neg[nNeg++] = i;
 109         }
 110         else
 111         {
 112             zero[nZero++] = i;
 113         }
 114     }
 115
 116     if (nNeg == 0)
 117     {
 118         return;
 119     }
 120
 121     if (nPos == 0)
 122     {
 123         decompTets.append(tetra);
 124         return;
 125     }
 126
 127     // Tetrahedron is split by plane.  Determine how it is split and how to
 128     // decompose the negative-side portion into tetrahedra (6 cases).
 129     scalar w0, w1, invCDiff;
 130     vector intp[4];
 131
 132     if (nPos == 3)
 133     {
 134         // +++-
 135         for (i = 0; i < nPos; ++i)
 136         {
 137             invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
 138
 139             w0 = -C[neg[0]] * invCDiff;
 140             w1 = +C[pos[i]] * invCDiff;
 141
 142             tetra[pos[i]] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
 143         }
 144
 145         decompTets.append(tetra);
 146     }
 147     else
 148     if (nPos == 2)
 149     {
 150         if (nNeg == 2)
 151         {
 152             // ++--
 153             for (i = 0; i < nPos; ++i)
 154             {
 155                 invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
 156
 157                 w0 = -C[neg[0]] * invCDiff;
 158                 w1 = +C[pos[i]] * invCDiff;
 159
 160                 intp[i] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
 161             }
 162
 163             for (i = 0; i < nNeg; ++i)
 164             {
 165                 invCDiff = (1.0 / (C[pos[i]] - C[neg[1]]));
 166
 167                 w0 = -C[neg[1]] * invCDiff;
 168                 w1 = +C[pos[i]] * invCDiff;
 169
 170                 intp[i+2] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[1]]);
 171             }
 172
 173             tetra[pos[0]] = intp[2];
 174             tetra[pos[1]] = intp[1];
 175
 176             decompTets.append(tetra);
 177
 178             tmpTetra[0] = tetra[neg[1]];
 179             tmpTetra[1] = intp[3];
 180             tmpTetra[2] = intp[2];
 181             tmpTetra[3] = intp[1];
 182
 183             decompTets.append(tmpTetra);
 184
 185             tmpTetra[0] = tetra[neg[0]];
 186             tmpTetra[1] = intp[0];
 187             tmpTetra[2] = intp[1];
 188             tmpTetra[3] = intp[2];
 189
 190             decompTets.append(tmpTetra);
 191         }
 192         else
 193         {
 194             // ++-0
 195             for (i = 0; i < nPos; ++i)
 196             {
 197                 invCDiff = (1.0 / (C[pos[i]] - C[neg[0]]));
 198
 199                 w0 = -C[neg[0]] * invCDiff;
 200                 w1 = +C[pos[i]] * invCDiff;
 201
 202                 tetra[pos[i]] = (w0 * tetra[pos[i]]) + (w1 * tetra[neg[0]]);
 203             }
 204
 205             decompTets.append(tetra);
 206         }
 207     }
 208     else
 209     if (nPos == 1)
 210     {
 211         if (nNeg == 3)
 212         {
 213             // +---
 214             for (i = 0; i < nNeg; ++i)
 215             {
 216                 invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
 217
 218                 w0 = -C[neg[i]] * invCDiff;
 219                 w1 = +C[pos[0]] * invCDiff;
 220
 221                 intp[i] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[i]]);
 222             }
 223
 224             tetra[pos[0]] = intp[0];
 225
 226             decompTets.append(tetra);
 227
 228             tmpTetra[0] = intp[0];
 229             tmpTetra[1] = tetra[neg[1]];
 230             tmpTetra[2] = tetra[neg[2]];
 231             tmpTetra[3] = intp[1];
 232
 233             decompTets.append(tmpTetra);
 234
 235             tmpTetra[0] = tetra[neg[2]];
 236             tmpTetra[1] = intp[1];
 237             tmpTetra[2] = intp[2];
 238             tmpTetra[3] = intp[0];
 239
 240             decompTets.append(tmpTetra);
 241         }
 242         else
 243         if (nNeg == 2)
 244         {
 245             // +--0
 246             for (i = 0; i < nNeg; ++i)
 247             {
 248                 invCDiff = (1.0 / (C[pos[0]] - C[neg[i]]));
 249
 250                 w0 = -C[neg[i]] * invCDiff;
 251                 w1 = +C[pos[0]] * invCDiff;
 252
 253                 intp[i] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[i]]);
 254             }
 255
 256             tetra[pos[0]] = intp[0];
 257
 258             decompTets.append(tetra);
 259
 260             tmpTetra[0] = intp[1];
 261             tmpTetra[1] = tetra[zero[0]];
 262             tmpTetra[2] = tetra[neg[1]];
 263             tmpTetra[3] = intp[0];
 264
 265             decompTets.append(tmpTetra);
 266         }
 267         else
 268         {
 269             // +-00
 270             invCDiff = (1.0 / (C[pos[0]] - C[neg[0]]));
 271
 272             w0 = -C[neg[0]] * invCDiff;
 273             w1 = +C[pos[0]] * invCDiff;
 274
 275             tetra[pos[0]] = (w0 * tetra[pos[0]]) + (w1 * tetra[neg[0]]);
 276
 277             decompTets.append(tetra);
 278         }
 279     }
 280 }
 281
 282
 283 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 284
 285 inline tetIntersection::tetIntersection(const FixedList<point, 4>& clipTet)
 286 :
 287     clipTet_(clipTet),
 288     clipTetMag_(0.0),
 289     inside_(10),
 290     allTets_(10)
 291 {
 292     // Pre-compute clipping planes
 293     computeClipPlanes();
 294 }
 295
 296
 297 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 298
 299 inline tetIntersection::~tetIntersection()
 300 {}
 301
 302
 303 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 304
 305 // Return magnitude of clipping tetrahedron
 306 inline scalar tetIntersection::clipTetMag() const
 307 {
 308     return clipTetMag_;
 309 }
 310
 311
 312 // Evaluate for intersections
 313 inline bool tetIntersection::evaluate(const FixedList<point, 4>& subjectTet)
 314 {
 315     // Clear lists
 316     inside_.clear();
 317     allTets_.clear();
 318
 319     // Add initial tetrahedron to list
 320     allTets_.append(subjectTet);
 321
 322     // Clip second tetrahedron against planes of clipping tetrahedron
 323     for (label i = 0; i < 4; i++)
 324     {
 325         forAll(allTets_, tetI)
 326         {
 327             splitAndDecompose(i, allTets_[tetI], inside_);
 328         }
 329
 330         // Prep for next clipping plane
 331         allTets_ = inside_;
 332         inside_.clear();
 333     }
 334
 335     return (allTets_.size() > 0);
 336 }
 337
 338
 339 // Return intersections
 340 inline const DynamicList<FixedList<point, 4> >&
 341 tetIntersection::getIntersection() const
 342 {
 343     return allTets_;
 344 }
 345
 346
 347 //- Evaluate and return volume / centroid
 348 inline void tetIntersection::getVolumeAndCentre
 349 (
 350     scalar& volume,
 351     vector& centre
 352 ) const
 353 {
 354     volume = 0.0;
 355     centre = vector::zero;
 356
 357     forAll(allTets_, tetI)
 358     {
 359         const FixedList<point, 4>& t = allTets_[tetI];
 360
 361         // Calculate volume (no check for orientation)
 362         scalar tV =
 363         (
 364             Foam::mag
 365             (
 366                 (1.0/6.0) *
 367                 (
 368                     ((t[1] - t[0]) ^ (t[2] - t[0])) & (t[3] - t[0])
 369                 )
 370             )
 371         );
 372
 373         // Calculate centroid
 374         vector tC = (0.25 * (t[0] + t[1] + t[2] + t[3]));
 375
 376         volume += tV;
 377         centre += (tV * tC);
 378     }
 379
 380     centre /= volume + VSMALL;
 381 }
 382
 383
 384 }
 385
 386 // ************************************************************************* //