src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/convexSetAlgorithm/tetIntersectionI.H

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