src/dynamicMesh/dynamicTopoFvMesh/convexSetAlgorithm/faceSetAlgorithm.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     faceSetAlgorithm
  26
  27 Description
  28     Class for convexSetAlgorithms pertaining to faces
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "faceSetAlgorithm.H"
  38
  39 #include "meshOps.H"
  40 #include "triFace.H"
  41 #include "polyMesh.H"
  42 #include "triIntersection.H"
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  50
  51 // Construct from components
  52 faceSetAlgorithm::faceSetAlgorithm
  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 faceSetAlgorithm::computeNormFactor(const label index) const
  79 {
  80     // Clear existing fields
  81     parents_.clear();
  82
  83     centres_.clear();
  84     weights_.clear();
  85
  86     // Compute refNorm and normFactor
  87     refNorm_ = newFaces_[index].normal(newPoints_);
  88     normFactor_ = mag(refNorm_);
  89
  90     // Normalize for later use
  91     refNorm_ /= normFactor_ + VSMALL;
  92
  93     // Compute a bounding box around the face
  94     box_ = boundBox(newFaces_[index].points(newPoints_), false);
  95
  96     vector minToXb = (box_.min() - box_.midpoint());
  97     vector maxToXb = (box_.max() - box_.midpoint());
  98
  99     // Scale it by a bit
 100     box_.min() += (1.5 * minToXb);
 101     box_.max() += (1.5 * maxToXb);
 102 }
 103
 104
 105 // Check whether the bounding box contains the entity
 106 bool faceSetAlgorithm::contains(const label index) const
 107 {
 108     // Fetch old face
 109     const face& checkFace = mesh_.faces()[index];
 110
 111     // Check if the bounding box contains any of the supplied points
 112     forAll(checkFace, pointI)
 113     {
 114         if (box_.contains(mesh_.points()[checkFace[pointI]]))
 115         {
 116             return true;
 117         }
 118     }
 119
 120     return false;
 121 }
 122
 123
 124 // Compute intersection
 125 bool faceSetAlgorithm::computeIntersection
 126 (
 127     const label newIndex,
 128     const label oldIndex,
 129     const label offset,
 130     bool output
 131 ) const
 132 {
 133     bool intersects = false;
 134
 135     const pointField& newPoints = newPoints_;
 136     const pointField& oldPoints = mesh_.points();
 137
 138     const face& newFace = newFaces_[newIndex];
 139     const face& oldFace = mesh_.faces()[oldIndex];
 140
 141     // Compute the normal for the old face
 142     vector oldNorm = oldFace.normal(oldPoints);
 143
 144     // Normalize
 145     oldNorm /= mag(oldNorm) + VSMALL;
 146
 147     if ((oldNorm & refNorm_) < 0.0)
 148     {
 149         // Opposite face orientation. Skip it.
 150         return false;
 151     }
 152
 153     // Check if decomposition is necessary
 154     if (oldFace.size() > 3 || newFace.size() > 3)
 155     {
 156         // Decompose new / old faces
 157         DynamicList<FixedList<point, 3> > clippingTris(15);
 158         DynamicList<FixedList<point, 3> > subjectTris(15);
 159
 160         label ntOld = 0, ntNew = 0;
 161
 162         if (oldFace.size() == 3)
 163         {
 164             // Add a new entry
 165             subjectTris.append(FixedList<point, 3>(vector::zero));
 166
 167             subjectTris[ntOld][0] = oldPoints[oldFace[0]];
 168             subjectTris[ntOld][1] = oldPoints[oldFace[1]];
 169             subjectTris[ntOld][2] = oldPoints[oldFace[2]];
 170
 171             ntOld++;
 172         }
 173         else
 174         {
 175             // Configure tris from oldFace
 176             vector ofCentre = oldFace.centre(oldPoints);
 177
 178             forAll(oldFace, pI)
 179             {
 180                 // Add a new entry
 181                 subjectTris.append(FixedList<point, 3>(vector::zero));
 182
 183                 subjectTris[ntOld][0] = oldPoints[oldFace[pI]];
 184                 subjectTris[ntOld][1] = oldPoints[oldFace.nextLabel(pI)];
 185                 subjectTris[ntOld][2] = ofCentre;
 186
 187                 ntOld++;
 188             }
 189         }
 190
 191         if (newFace.size() == 3)
 192         {
 193             // Add a new entry
 194             clippingTris.append(FixedList<point, 3>(vector::zero));
 195
 196             clippingTris[ntNew][0] = newPoints[newFace[0]];
 197             clippingTris[ntNew][1] = newPoints[newFace[1]];
 198             clippingTris[ntNew][2] = newPoints[newFace[2]];
 199
 200             ntNew++;
 201         }
 202         else
 203         {
 204             // Configure tris from newFace
 205             vector nfCentre = newFace.centre(newPoints);
 206
 207             forAll(newFace, pI)
 208             {
 209                 // Add a new entry
 210                 clippingTris.append(FixedList<point, 3>(vector::zero));
 211
 212                 clippingTris[ntNew][0] = newPoints[newFace[pI]];
 213                 clippingTris[ntNew][1] = newPoints[newFace.nextLabel(pI)];
 214                 clippingTris[ntNew][2] = nfCentre;
 215
 216                 ntNew++;
 217             }
 218         }
 219
 220         // Accumulate area / centroid over all intersections
 221         bool foundIntersect = false;
 222
 223         scalar totalArea = 0.0;
 224         vector totalCentre = vector::zero;
 225
 226         // Loop through all clipping tris
 227         forAll(clippingTris, i)
 228         {
 229             // Initialize the intersector
 230             triIntersection intersector(clippingTris[i]);
 231
 232             // Test for intersection and evaluate
 233             // against all subject tris
 234             forAll(subjectTris, j)
 235             {
 236                 intersects = intersector.evaluate(subjectTris[j]);
 237
 238                 if (intersects)
 239                 {
 240                     scalar area = 0.0;
 241                     vector centre = vector::zero;
 242
 243                     // Fetch area and centre
 244                     intersector.getAreaAndCentre(area, centre);
 245
 246                     // Accumulate area / centroid
 247                     totalArea += area;
 248                     totalCentre += (area * centre);
 249
 250                     foundIntersect = true;
 251
 252                     if (output)
 253                     {
 254                         writeVTK
 255                         (
 256                             "polygonIntersect_"
 257                           + Foam::name(newIndex)
 258                           + '_'
 259                           + Foam::name(oldIndex)
 260                           + '_' + Foam::name(i)
 261                           + '_' + Foam::name(j),
 262                             intersector.getIntersection()
 263                         );
 264                     }
 265                 }
 266             }
 267         }
 268
 269         // Size-up the internal lists
 270         if (foundIntersect && !output)
 271         {
 272             // Normalize centre
 273             totalCentre /= totalArea + VSMALL;
 274
 275             // Normalize and check if this is worth it
 276             if (mag(totalArea/normFactor_) > SMALL)
 277             {
 278                 // Size-up the internal lists
 279                 meshOps::sizeUpList((oldIndex - offset), parents_);
 280                 meshOps::sizeUpList(totalArea, weights_);
 281                 meshOps::sizeUpList(totalCentre, centres_);
 282
 283                 intersects = true;
 284             }
 285             else
 286             {
 287                 intersects = false;
 288             }
 289         }
 290     }
 291     else
 292     {
 293         // Configure points for clipping triangle
 294         FixedList<point, 3> clippingTri(vector::zero);
 295
 296         // Fill in points
 297         clippingTri[0] = newPoints[newFace[0]];
 298         clippingTri[1] = newPoints[newFace[1]];
 299         clippingTri[2] = newPoints[newFace[2]];
 300
 301         // Configure points for subject triangle
 302         FixedList<point, 3> subjectTri(vector::zero);
 303
 304         // Fill in points
 305         subjectTri[0] = oldPoints[oldFace[0]];
 306         subjectTri[1] = oldPoints[oldFace[1]];
 307         subjectTri[2] = oldPoints[oldFace[2]];
 308
 309         // Initialize the intersector
 310         triIntersection intersector(clippingTri);
 311
 312         // Test for intersection and evaluate
 313         intersects = intersector.evaluate(subjectTri);
 314
 315         if (intersects)
 316         {
 317             scalar area;
 318             vector centre;
 319
 320             // Fetch area and centre
 321             intersector.getAreaAndCentre(area, centre);
 322
 323             // Normalize and check if this is worth it
 324             if (mag(area/normFactor_) > SMALL)
 325             {
 326                 if (output)
 327                 {
 328                     writeVTK
 329                     (
 330                         "triIntersectNew_"
 331                       + Foam::name(newIndex),
 332                         List<FixedList<point, 3> >(1, clippingTri)
 333                     );
 334
 335                     writeVTK
 336                     (
 337                         "triIntersectOld_"
 338                       + Foam::name(newIndex)
 339                       + '_'
 340                       + Foam::name(oldIndex),
 341                         List<FixedList<point, 3> >(1, subjectTri)
 342                     );
 343
 344                     writeVTK
 345                     (
 346                         "triIntersect_"
 347                       + Foam::name(newIndex)
 348                       + '_'
 349                       + Foam::name(oldIndex),
 350                         intersector.getIntersection()
 351                     );
 352                 }
 353
 354                 // Size-up the internal lists
 355                 meshOps::sizeUpList((oldIndex - offset), parents_);
 356                 meshOps::sizeUpList(area, weights_);
 357                 meshOps::sizeUpList(centre, centres_);
 358             }
 359             else
 360             {
 361                 intersects = false;
 362             }
 363         }
 364     }
 365
 366     return intersects;
 367 }
 368
 369
 370 //- Write out tris as a VTK
 371 void faceSetAlgorithm::writeVTK
 372 (
 373     const word& name,
 374     const List<FixedList<point, 3> >& triList
 375 ) const
 376 {
 377     // Fill up all points
 378     label pI = 0;
 379
 380     List<face> allTris(triList.size(), face(3));
 381     pointField allPoints(3 * triList.size());
 382
 383     forAll(triList, triI)
 384     {
 385         allTris[triI][0] = pI;
 386         allPoints[pI++] = triList[triI][0];
 387
 388         allTris[triI][1] = pI;
 389         allPoints[pI++] = triList[triI][1];
 390
 391         allTris[triI][2] = pI;
 392         allPoints[pI++] = triList[triI][2];
 393     }
 394
 395     // Write out in face-to-node addressing
 396     meshOps::writeVTK
 397     (
 398         mesh_,
 399         name,
 400         identity(triList.size()),
 401         2,
 402         allPoints,
 403         List<edge>(0),
 404         allTris,
 405         List<cell>(0),
 406         List<label>(0)
 407     );
 408 }
 409
 410
 411 } // End namespace Foam
 412
 413 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //