src/dynamicMesh/dynamicTopoFvMesh/convexSetAlgorithm/cellSetAlgorithm.C

   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 Class
  25     cellSetAlgorithm
  26
  27 Description
  28     Class for convexSetAlgorithms pertaining to cells
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "cellSetAlgorithm.H"
  38
  39 #include "meshOps.H"
  40 #include "triFace.H"
  41 #include "polyMesh.H"
  42 #include "tetIntersection.H"
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  50
  51 // Construct from components
  52 cellSetAlgorithm::cellSetAlgorithm
  53 (
  54     const polyMesh& mesh,
  55     const pointField& newPoints,
  56     const UList<edge>& newEdges,
  57     const UList<face>& newFaces,
  58     const UList<cell>& newCells,
  59     const UList<label>& newOwner,
  60     const UList<label>& newNeighbour
  61 )
  62 :
  63     convexSetAlgorithm
  64     (
  65         mesh,
  66         newPoints,
  67         newEdges,
  68         newFaces,
  69         newCells,
  70         newOwner,
  71         newNeighbour
  72     )
  73 {}
  74
  75
  76 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  77
  78 void cellSetAlgorithm::computeNormFactor(const label index) const
  79 {
  80     // Clear existing fields
  81     parents_.clear();
  82
  83     centres_.clear();
  84     weights_.clear();
  85
  86     // Compute volume / centre (using refNorm_ as centre)
  87     meshOps::cellCentreAndVolume
  88     (
  89         index,
  90         newPoints_,
  91         newFaces_,
  92         newCells_,
  93         newOwner_,
  94         refNorm_,
  95         normFactor_
  96     );
  97
  98     // Compute a bounding box around the cell
  99     box_ = boundBox(newCells_[index].points(newFaces_, newPoints_), false);
 100
 101     vector minToXb = (box_.min() - box_.midpoint());
 102     vector maxToXb = (box_.max() - box_.midpoint());
 103
 104     // Scale it by a bit
 105     box_.min() += (1.5 * minToXb);
 106     box_.max() += (1.5 * maxToXb);
 107 }
 108
 109
 110 // Check whether the bounding box contains the entity
 111 bool cellSetAlgorithm::contains(const label index) const
 112 {
 113     // Fetch old points
 114     const labelListList& cellPoints = mesh_.cellPoints();
 115     const labelList& checkCell = cellPoints[index];
 116
 117     // Check if the bounding box contains any of the supplied points
 118     forAll(checkCell, pointI)
 119     {
 120         if (box_.contains(mesh_.points()[checkCell[pointI]]))
 121         {
 122             return true;
 123         }
 124     }
 125
 126     return false;
 127 }
 128
 129
 130 // Compute intersection
 131 bool cellSetAlgorithm::computeIntersection
 132 (
 133     const label newIndex,
 134     const label oldIndex,
 135     const label offset,
 136     bool output
 137 ) const
 138 {
 139     bool intersects = false;
 140
 141     const pointField& newPoints = newPoints_;
 142     const pointField& oldPoints = mesh_.points();
 143
 144     const cell& newCell = newCells_[newIndex];
 145     const cell& oldCell = mesh_.cells()[oldIndex];
 146
 147     // Check if decomposition is necessary
 148     if (oldCell.size() > 4 || newCell.size() > 4)
 149     {
 150         // Decompose new / old cells
 151         DynamicList<FixedList<point, 4> > clippingTets(15);
 152         DynamicList<FixedList<point, 4> > subjectTets(15);
 153
 154         label ntOld = 0, ntNew = 0;
 155         vector oldCentre = vector::zero, newCentre = vector::zero;
 156
 157         // Configure tets from oldCell
 158         forAll(oldCell, faceI)
 159         {
 160             const face& oldFace = mesh_.faces()[oldCell[faceI]];
 161
 162             vector fCentre = oldFace.centre(oldPoints);
 163
 164             if (oldFace.size() == 3)
 165             {
 166                 // Add a new entry
 167                 subjectTets.append(FixedList<point, 4>(vector::zero));
 168
 169                 subjectTets[ntOld][0] = oldPoints[oldFace[0]];
 170                 subjectTets[ntOld][1] = oldPoints[oldFace[1]];
 171                 subjectTets[ntOld][2] = oldPoints[oldFace[2]];
 172
 173                 ntOld++;
 174             }
 175             else
 176             {
 177                 forAll(oldFace, pI)
 178                 {
 179                     // Add a new entry
 180                     subjectTets.append(FixedList<point, 4>(vector::zero));
 181
 182                     subjectTets[ntOld][0] = oldPoints[oldFace[pI]];
 183                     subjectTets[ntOld][1] = oldPoints[oldFace.nextLabel(pI)];
 184                     subjectTets[ntOld][2] = fCentre;
 185
 186                     ntOld++;
 187                 }
 188             }
 189
 190             oldCentre += fCentre;
 191         }
 192
 193         // Configure tets from newCell
 194         forAll(newCell, faceI)
 195         {
 196             const face& newFace = newFaces_[newCell[faceI]];
 197
 198             vector fCentre = newFace.centre(newPoints);
 199
 200             if (newFace.size() == 3)
 201             {
 202                 // Add a new entry
 203                 clippingTets.append(FixedList<point, 4>(vector::zero));
 204
 205                 clippingTets[ntNew][0] = newPoints[newFace[0]];
 206                 clippingTets[ntNew][1] = newPoints[newFace[1]];
 207                 clippingTets[ntNew][2] = newPoints[newFace[2]];
 208
 209                 ntNew++;
 210             }
 211             else
 212             {
 213                 forAll(newFace, pI)
 214                 {
 215                     // Add a new entry
 216                     clippingTets.append(FixedList<point, 4>(vector::zero));
 217
 218                     clippingTets[ntNew][0] = newPoints[newFace[pI]];
 219                     clippingTets[ntNew][1] = newPoints[newFace.nextLabel(pI)];
 220                     clippingTets[ntNew][2] = fCentre;
 221
 222                     ntNew++;
 223                 }
 224             }
 225
 226             newCentre += fCentre;
 227         }
 228
 229         oldCentre /= oldCell.size();
 230         newCentre /= newCell.size();
 231
 232         // Fill-in last points for all tets
 233         forAll(subjectTets, i)
 234         {
 235             subjectTets[i][3] = oldCentre;
 236         }
 237
 238         forAll(clippingTets, i)
 239         {
 240             clippingTets[i][3] = newCentre;
 241         }
 242
 243         // Accumulate volume / centroid over all intersections
 244         bool foundIntersect = false;
 245
 246         scalar totalVolume = 0.0;
 247         vector totalCentre = vector::zero;
 248
 249         // Loop through all clipping tets
 250         forAll(clippingTets, i)
 251         {
 252             // Initialize the intersector
 253             tetIntersection intersector(clippingTets[i]);
 254
 255             // Test for intersection and evaluate
 256             // against all subject tets
 257             forAll(subjectTets, j)
 258             {
 259                 intersects = intersector.evaluate(subjectTets[j]);
 260
 261                 if (intersects)
 262                 {
 263                     scalar volume = 0.0;
 264                     vector centre = vector::zero;
 265
 266                     // Fetch volume and centre
 267                     intersector.getVolumeAndCentre(volume, centre);
 268
 269                     // Accumulate volume / centroid
 270                     totalVolume += volume;
 271                     totalCentre += (volume * centre);
 272
 273                     foundIntersect = true;
 274
 275                     if (output)
 276                     {
 277                         writeVTK
 278                         (
 279                             "polyhedronIntersect_"
 280                           + Foam::name(newIndex)
 281                           + '_'
 282                           + Foam::name(oldIndex)
 283                           + '_' + Foam::name(i)
 284                           + '_' + Foam::name(j),
 285                             intersector.getIntersection()
 286                         );
 287                     }
 288                 }
 289             }
 290         }
 291
 292         // Size-up the internal lists
 293         if (foundIntersect && !output)
 294         {
 295             // Normalize centre
 296             totalCentre /= totalVolume + VSMALL;
 297
 298             // Normalize and check if this is worth it
 299             if (mag(totalVolume/normFactor_) > SMALL)
 300             {
 301                 // Size-up the internal lists
 302                 meshOps::sizeUpList((oldIndex - offset), parents_);
 303                 meshOps::sizeUpList(totalVolume, weights_);
 304                 meshOps::sizeUpList(totalCentre, centres_);
 305
 306                 intersects = true;
 307             }
 308             else
 309             {
 310                 intersects = false;
 311             }
 312         }
 313     }
 314     else
 315     {
 316         // Configure points for clipping tetrahedron
 317         FixedList<point, 4> clippingTet(vector::zero);
 318
 319         // Configure the clipping tetrahedron.
 320         const face& firstNewFace = newFaces_[newCell[0]];
 321         const face& secondNewFace = newFaces_[newCell[1]];
 322
 323         // Find the isolated point
 324         label fourthNewPoint =
 325         (
 326             meshOps::findIsolatedPoint
 327             (
 328                 firstNewFace,
 329                 secondNewFace
 330             )
 331         );
 332
 333         // Fill in points
 334         clippingTet[0] = newPoints[firstNewFace[0]];
 335         clippingTet[1] = newPoints[firstNewFace[1]];
 336         clippingTet[2] = newPoints[firstNewFace[2]];
 337         clippingTet[3] = newPoints[fourthNewPoint];
 338
 339         // Configure points for subject tetrahedron
 340         FixedList<point, 4> subjectTet(vector::zero);
 341
 342         // Configure the subject tetrahedron.
 343         const face& firstOldFace = mesh_.faces()[oldCell[0]];
 344         const face& secondOldFace = mesh_.faces()[oldCell[1]];
 345
 346         // Find the isolated point
 347         label fourthOldPoint =
 348         (
 349             meshOps::findIsolatedPoint
 350             (
 351                 firstOldFace,
 352                 secondOldFace
 353             )
 354         );
 355
 356         // Fill in points
 357         subjectTet[0] = oldPoints[firstOldFace[0]];
 358         subjectTet[1] = oldPoints[firstOldFace[1]];
 359         subjectTet[2] = oldPoints[firstOldFace[2]];
 360         subjectTet[3] = oldPoints[fourthOldPoint];
 361
 362         // Initialize the intersector
 363         tetIntersection intersector(clippingTet);
 364
 365         // Test for intersection and evaluate
 366         intersects = intersector.evaluate(subjectTet);
 367
 368         if (intersects)
 369         {
 370             scalar volume = 0.0;
 371             vector centre = vector::zero;
 372
 373             // Fetch volume and centre
 374             intersector.getVolumeAndCentre(volume, centre);
 375
 376             // Normalize and check if this is worth it
 377             if (mag(volume/normFactor_) > SMALL)
 378             {
 379                 if (output)
 380                 {
 381                     writeVTK
 382                     (
 383                         "tetIntersectNew_"
 384                       + Foam::name(newIndex),
 385                         List<FixedList<point, 4> >(1, clippingTet)
 386                     );
 387
 388                     writeVTK
 389                     (
 390                         "tetIntersectOld_"
 391                       + Foam::name(newIndex)
 392                       + '_'
 393                       + Foam::name(oldIndex),
 394                         List<FixedList<point, 4> >(1, subjectTet)
 395                     );
 396
 397                     writeVTK
 398                     (
 399                         "tetIntersect_"
 400                       + Foam::name(newIndex)
 401                       + '_'
 402                       + Foam::name(oldIndex),
 403                         intersector.getIntersection()
 404                     );
 405                 }
 406
 407                 // Size-up the internal lists
 408                 meshOps::sizeUpList((oldIndex - offset), parents_);
 409                 meshOps::sizeUpList(volume, weights_);
 410                 meshOps::sizeUpList(centre, centres_);
 411             }
 412             else
 413             {
 414                 intersects = false;
 415             }
 416         }
 417     }
 418
 419     return intersects;
 420 }
 421
 422
 423 //- Write out tets as a VTK
 424 void cellSetAlgorithm::writeVTK
 425 (
 426     const word& name,
 427     const List<FixedList<point, 4> >& tetList
 428 ) const
 429 {
 430     // Fill up all points
 431     label pI = 0;
 432
 433     List<cell> allTets(tetList.size(), cell(4));
 434     pointField allPoints(4 * tetList.size());
 435
 436     forAll(tetList, tetI)
 437     {
 438         allTets[tetI][0] = pI;
 439         allPoints[pI++] = tetList[tetI][0];
 440
 441         allTets[tetI][1] = pI;
 442         allPoints[pI++] = tetList[tetI][1];
 443
 444         allTets[tetI][2] = pI;
 445         allPoints[pI++] = tetList[tetI][2];
 446
 447         allTets[tetI][3] = pI;
 448         allPoints[pI++] = tetList[tetI][3];
 449     }
 450
 451     // Write out in cell-to-node addressing
 452     meshOps::writeVTK
 453     (
 454         mesh_,
 455         name,
 456         identity(tetList.size()),
 457         3,
 458         allPoints,
 459         List<edge>(0),
 460         List<face>(0),
 461         allTets,
 462         List<label>(0)
 463     );
 464 }
 465
 466
 467 } // End namespace Foam
 468
 469 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //