src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.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 Class
  26     meshOps
  27
  28 Description
  29     Various utility functions that perform mesh-related operations.
  30
  31 Author
  32     Sandeep Menon
  33     University of Massachusetts Amherst
  34     All rights reserved
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "Pair.H"
  39 #include "meshOps.H"
  40 #include "ListOps.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 namespace meshOps
  48 {
  49
  50 // Utility method to find the common edge between two faces.
  51 inline bool findCommonEdge
  52 (
  53     const label first,
  54     const label second,
  55     const UList<labelList>& faceEdges,
  56     label& common
  57 )
  58 {
  59     bool found = false;
  60
  61     const labelList& fEi = faceEdges[first];
  62     const labelList& fEj = faceEdges[second];
  63
  64     forAll(fEi, edgeI)
  65     {
  66         forAll(fEj, edgeJ)
  67         {
  68             if (fEi[edgeI] == fEj[edgeJ])
  69             {
  70                 common = fEi[edgeI];
  71
  72                 found = true;
  73                 break;
  74             }
  75         }
  76
  77         if (found)
  78         {
  79             break;
  80         }
  81     }
  82
  83     return found;
  84 }
  85
  86
  87 // Utility method to find the interior (quad) / boundary (tri) faces
  88 // for an input quad-face and adjacent triangle-prism cell.
  89 inline void findPrismFaces
  90 (
  91     const label fIndex,
  92     const label cIndex,
  93     const UList<face>& faces,
  94     const UList<cell>& cells,
  95     const UList<label>& neighbour,
  96     FixedList<face,2>& bdyf,
  97     FixedList<label,2>& bidx,
  98     FixedList<face,2>& intf,
  99     FixedList<label,2>& iidx
 100 )
 101 {
 102     label indexO = 0, indexI = 0;
 103
 104     const cell& c = cells[cIndex];
 105
 106     forAll(c, i)
 107     {
 108         label faceIndex = c[i];
 109
 110         // Don't count the face under consideration
 111         if (faceIndex != fIndex)
 112         {
 113             const face& fi = faces[faceIndex];
 114
 115             if (neighbour[faceIndex] == -1)
 116             {
 117                 if (fi.size() == 3)
 118                 {
 119                     // Triangular face on the boundary
 120                     bidx[indexO] = faceIndex;
 121                     bdyf[indexO++] = fi;
 122                 }
 123                 else
 124                 {
 125                     // This seems to be a non-triangular face on the boundary
 126                     // Consider this as "interior" and move on
 127                     iidx[indexI] = faceIndex;
 128                     intf[indexI++] = fi;
 129                 }
 130             }
 131             else
 132             {
 133                 // Face on the interior
 134                 iidx[indexI] = faceIndex;
 135                 intf[indexI++] = fi;
 136             }
 137         }
 138     }
 139 }
 140
 141
 142 // Utility method to find the isolated point given two triangular faces.
 143 //  - Returns the point on checkFace that does not belong to baseFace.
 144 inline label findIsolatedPoint
 145 (
 146     const face& baseFace,
 147     const face& checkFace
 148 )
 149 {
 150     // Get the fourth point
 151     forAll(checkFace, pointI)
 152     {
 153         if
 154         (
 155             checkFace[pointI] != baseFace[0] &&
 156             checkFace[pointI] != baseFace[1] &&
 157             checkFace[pointI] != baseFace[2]
 158         )
 159         {
 160             return checkFace[pointI];
 161         }
 162     }
 163
 164     return -1;
 165 }
 166
 167
 168 // Utility method to find the isolated point on a triangular face
 169 // that doesn't lie on the specified edge. Also returns the point next to it.
 170 inline void findIsolatedPoint
 171 (
 172     const face& f,
 173     const edge& e,
 174     label& ptIndex,
 175     label& nextPtIndex
 176 )
 177 {
 178     // Check the first point
 179     if ( f[0] != e.start() && f[0] != e.end() )
 180     {
 181         ptIndex = f[0];
 182         nextPtIndex = f[1];
 183         return;
 184     }
 185
 186     // Check the second point
 187     if ( f[1] != e.start() && f[1] != e.end() )
 188     {
 189         ptIndex = f[1];
 190         nextPtIndex = f[2];
 191         return;
 192     }
 193
 194     // Check the third point
 195     if ( f[2] != e.start() && f[2] != e.end() )
 196     {
 197         ptIndex = f[2];
 198         nextPtIndex = f[0];
 199         return;
 200     }
 201
 202     // This bit should never happen.
 203     FatalErrorIn
 204     (
 205         "inline void meshOps::findIsolatedPoint"
 206         "(const face&, const edge&, label&, label&)"
 207     )
 208         << "Cannot find isolated point in face " << f << endl
 209         << " Using edge: " << e
 210         << abort(FatalError);
 211 }
 212
 213
 214 // Given a face and cell index, find the apex point for a tet cell.
 215 inline label tetApexPoint
 216 (
 217     const label cIndex,
 218     const label fIndex,
 219     const UList<face>& faces,
 220     const UList<cell>& cells
 221 )
 222 {
 223     label apexPoint = -1;
 224     bool foundApex = false;
 225
 226     const cell& cellToCheck = cells[cIndex];
 227     const face& baseFace = faces[fIndex];
 228
 229     forAll(cellToCheck, faceI)
 230     {
 231         const face& faceToCheck = faces[cellToCheck[faceI]];
 232
 233         apexPoint = findIsolatedPoint(baseFace, faceToCheck);
 234
 235         if (apexPoint > -1)
 236         {
 237             foundApex = true;
 238             break;
 239         }
 240     }
 241
 242     if (!foundApex)
 243     {
 244         Pout<< "fIndex: " << fIndex << ":: " << faces[fIndex] << endl;
 245         Pout<< "cIndex: " << cIndex << ":: " << cellToCheck << endl;
 246
 247         forAll(cellToCheck, faceI)
 248         {
 249             Pout<< '\t' << cellToCheck[faceI] << ":: "
 250                 << faces[cellToCheck[faceI]] << endl;
 251         }
 252
 253         FatalErrorIn
 254         (
 255             "\n\n"
 256             "inline label meshOps::tetApexPoint\n"
 257             "(\n"
 258             "    const label cIndex,\n"
 259             "    const label fIndex,\n"
 260             "    const UList<face>& faces,\n"
 261             "    const UList<cell>& cells\n"
 262             ")\n"
 263         )
 264             << "Could not find an apex point in cell: " << cIndex
 265             << " given face: " << fIndex
 266             << abort(FatalError);
 267     }
 268
 269     return apexPoint;
 270 }
 271
 272
 273 // Given a cell index, find the centroid / volume of a cell.
 274 //  - If checking is enabled, return whether the cell is closed
 275 inline void cellCentreAndVolume
 276 (
 277     const label cIndex,
 278     const vectorField& points,
 279     const UList<face>& faces,
 280     const UList<cell>& cells,
 281     const UList<label>& owner,
 282     vector& centre,
 283     scalar& volume
 284 )
 285 {
 286     // Reset inputs
 287     volume = 0.0;
 288     centre = vector::zero;
 289
 290     const cell& cellToCheck = cells[cIndex];
 291
 292     // Average face-centres to get an estimate centroid
 293     vector cEst(vector::zero), fCentre(vector::zero);
 294
 295     forAll(cellToCheck, faceI)
 296     {
 297         const face& checkFace = faces[cellToCheck[faceI]];
 298
 299         cEst += checkFace.centre(points);
 300     }
 301
 302     cEst /= cellToCheck.size();
 303
 304     forAll(cellToCheck, faceI)
 305     {
 306         const face& checkFace = faces[cellToCheck[faceI]];
 307
 308         if (checkFace.size() == 3)
 309         {
 310             tetPointRef tpr
 311             (
 312                 points[checkFace[2]],
 313                 points[checkFace[1]],
 314                 points[checkFace[0]],
 315                 cEst
 316             );
 317
 318             scalar tetVol = tpr.mag();
 319
 320             // Flip if necessary
 321             if (owner[cellToCheck[faceI]] != cIndex)
 322             {
 323                 tetVol *= -1.0;
 324             }
 325
 326             // Accumulate volume-weighted tet centre
 327             centre += tetVol*tpr.centre();
 328
 329             // Accumulate tet volume
 330             volume += tetVol;
 331         }
 332         else
 333         {
 334             forAll(checkFace, pI)
 335             {
 336                 tetPointRef tpr
 337                 (
 338                     points[checkFace[pI]],
 339                     points[checkFace.prevLabel(pI)],
 340                     checkFace.centre(points),
 341                     cEst
 342                 );
 343
 344                 scalar tetVol = tpr.mag();
 345
 346                 // Flip if necessary
 347                 if (owner[cellToCheck[faceI]] != cIndex)
 348                 {
 349                     tetVol *= -1.0;
 350                 }
 351
 352                 // Accumulate volume-weighted tet centre
 353                 centre += tetVol*tpr.centre();
 354
 355                 // Accumulate tet volume
 356                 volume += tetVol;
 357             }
 358         }
 359     }
 360
 361     centre /= volume + VSMALL;
 362 }
 363
 364
 365 // Determine whether a segment intersects a triangular face
 366 inline bool segmentTriFaceIntersection
 367 (
 368     const triPointRef& faceToCheck,
 369     const linePointRef& edgeToCheck,
 370     vector& intPoint
 371 )
 372 {
 373     // Fetch references
 374     const vector& p1 = edgeToCheck.start();
 375     const vector& p2 = edgeToCheck.end();
 376
 377     // Define face-normal
 378     vector n = faceToCheck.normal();
 379     n /= mag(n) + VSMALL;
 380
 381     // Compute uValue
 382     scalar numerator = n & (faceToCheck.a() - p1);
 383     scalar denominator = n & (p2 - p1);
 384
 385     // Check if the edge is parallel to the face
 386     if (mag(denominator) < VSMALL)
 387     {
 388         return false;
 389     }
 390
 391     scalar u = (numerator / denominator);
 392
 393     // Check for intersection along line.
 394     if ((u >= 0.0) && (u <= 1.0))
 395     {
 396         // Compute point of intersection
 397         intPoint = p1 + u*(p2 - p1);
 398
 399         // Also make sure that intPoint lies within the face
 400         if (pointInTriFace(faceToCheck, intPoint, false))
 401         {
 402             return true;
 403         }
 404     }
 405
 406     // Failed to fall within edge-bounds, or within face
 407     return false;
 408 }
 409
 410
 411 // Determine whether the particular point lies
 412 // inside the given triangular face
 413 inline bool pointInTriFace
 414 (
 415     const triPointRef& faceToCheck,
 416     const vector& cP,
 417     bool testCoplanarity
 418 )
 419 {
 420     // Copy inputs
 421     const vector& a = faceToCheck.a();
 422     const vector& b = faceToCheck.b();
 423     const vector& c = faceToCheck.c();
 424
 425     // Compute the normal
 426     vector nf = faceToCheck.normal();
 427
 428     if ( ((0.5 * ((b - a)^(cP - a))) & nf) < 0.0)
 429     {
 430         return false;
 431     }
 432
 433     if ( ((0.5 * ((c - b)^(cP - b))) & nf) < 0.0)
 434     {
 435         return false;
 436     }
 437
 438     if ( ((0.5 * ((a - c)^(cP - c))) & nf) < 0.0)
 439     {
 440         return false;
 441     }
 442
 443     if (testCoplanarity)
 444     {
 445         vector ftp(a - cP);
 446
 447         // Normalize vectors
 448         nf /= mag(nf) + VSMALL;
 449         ftp /= mag(ftp) + VSMALL;
 450
 451         if (mag(ftp & nf) > VSMALL)
 452         {
 453             return false;
 454         }
 455     }
 456
 457     // Passed test with all edges
 458     return true;
 459 }
 460
 461
 462 // Parallel blocking send
 463 inline void pWrite
 464 (
 465     const label toID,
 466     const label& data
 467 )
 468 {
 469     OPstream::write
 470     (
 471         Pstream::blocking,
 472         toID,
 473         reinterpret_cast<const char*>(&data),
 474         sizeof(label)
 475     );
 476 }
 477
 478
 479 // Parallel blocking receive
 480 inline void pRead
 481 (
 482     const label fromID,
 483     label& data
 484 )
 485 {
 486     IPstream::read
 487     (
 488         Pstream::blocking,
 489         fromID,
 490         reinterpret_cast<char*>(&data),
 491         sizeof(label)
 492     );
 493 }
 494
 495
 496 inline void waitForBuffers()
 497 {
 498     if (Pstream::parRun())
 499     {
 500         OPstream::waitRequests();
 501         IPstream::waitRequests();
 502     }
 503 }
 504
 505
 506 // Method to insert a label between two labels in a list
 507 // Assumes that all labels are unique.
 508 inline void insertLabel
 509 (
 510     const label newLabel,
 511     const label labelA,
 512     const label labelB,
 513     labelList& list
 514 )
 515 {
 516     // Create a new list
 517     bool found = false;
 518     label origSize = list.size();
 519     labelList newList(origSize + 1);
 520
 521     label index = 0, nextI = -1;
 522
 523     // Start a linear search
 524     forAll(list, itemI)
 525     {
 526         newList[index++] = list[itemI];
 527
 528         nextI = list.fcIndex(itemI);
 529
 530         if
 531         (
 532             (
 533                 (list[itemI] == labelA && list[nextI] == labelB) ||
 534                 (list[itemI] == labelB && list[nextI] == labelA)
 535             ) &&
 536            !found
 537         )
 538         {
 539             found = true;
 540             newList[index++] = newLabel;
 541         }
 542     }
 543
 544     if (!found)
 545     {
 546         FatalErrorIn
 547         (
 548             "inline void meshOps::insertLabel"
 549             "(const label, const label, const label, labelList&)"
 550         )   << nl << "Cannot insert " << newLabel
 551             << " in list: " << list << nl
 552             << " Labels: "
 553             << labelA << " and " << labelB
 554             << " were not found in sequence."
 555             << abort(FatalError);
 556     }
 557
 558     // Transfer the list
 559     list.transfer(newList);
 560 }
 561
 562
 563 // Utility method to replace a label in a given list
 564 inline void replaceLabel
 565 (
 566      const label original,
 567      const label replacement,
 568      labelList& list
 569 )
 570 {
 571     label index = -1;
 572
 573     if ((index = findIndex(list, original)) > -1)
 574     {
 575         list[index] = replacement;
 576     }
 577     else
 578     {
 579         FatalErrorIn
 580         (
 581             "inline void label meshOps::replaceLabel"
 582             "(const label, const label, labelList&)"
 583         )   << nl << "Cannot find " << original
 584             << " in list: " << list << nl
 585             << " Label: " << replacement
 586             << " was not used in replacement."
 587             << abort(FatalError);
 588     }
 589 }
 590
 591
 592 } // End namespace meshOps
 593
 594
 595 } // End namespace Foam
 596
 597 // ************************************************************************* //