src/dynamicMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.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 Class
  25     dynamicTopoFvMesh
  26
  27 Description
  28     An implementation of dynamic changes to mesh-topology
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "Stack.H"
  38 #include "meshOps.H"
  39 #include "tetPointRef.H"
  40 #include "linePointRef.H"
  41 #include "lengthScaleEstimator.H"
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  49
  50 // Return if mesh is 2D
  51 inline bool dynamicTopoFvMesh::is2D() const
  52 {
  53     return twoDMesh_;
  54 }
  55
  56
  57 // Return if mesh is 3D
  58 inline bool dynamicTopoFvMesh::is3D() const
  59 {
  60     return !twoDMesh_;
  61 }
  62
  63
  64 // Return a non-const reference to the lengthScaleEstimator
  65 inline lengthScaleEstimator&
  66 dynamicTopoFvMesh::lengthEstimator()
  67 {
  68     if (!lengthEstimator_.valid())
  69     {
  70         FatalErrorIn
  71         (
  72             "inline lengthScaleEstimator& "
  73             "dynamicTopoFvMesh::lengthEstimator()"
  74         ) << nl
  75           << " Invalid request. Length scale estimator was not loaded. "
  76           << abort(FatalError);
  77     }
  78
  79     return lengthEstimator_();
  80 }
  81
  82
  83 // Return a const reference to the lengthScaleEstimator
  84 inline const lengthScaleEstimator&
  85 dynamicTopoFvMesh::lengthEstimator() const
  86 {
  87     if (!lengthEstimator_.valid())
  88     {
  89         FatalErrorIn
  90         (
  91             "inline const lengthScaleEstimator& "
  92             "dynamicTopoFvMesh::lengthEstimator() const"
  93         ) << nl
  94           << " Invalid request. Length scale estimator was not loaded. "
  95           << abort(FatalError);
  96     }
  97
  98     return lengthEstimator_();
  99 }
 100
 101
 102 // Return a const reference to the multiThreader
 103 inline const multiThreader&
 104 dynamicTopoFvMesh::threader() const
 105 {
 106     if (!threader_.valid())
 107     {
 108         FatalErrorIn
 109         (
 110             "inline const multiThreader& "
 111             "dynamicTopoFvMesh::threader() const"
 112         ) << nl
 113           << " Invalid request. multiThreader was not loaded. "
 114           << abort(FatalError);
 115     }
 116
 117     return threader_();
 118 }
 119
 120
 121 // Return a reference to the entity mutexes.
 122 // The index 'entity' ranges from 0 to 3 for point/edge/face/cell.
 123 inline const Mutex& dynamicTopoFvMesh::entityMutex
 124 (
 125     const label entity
 126 ) const
 127 {
 128     return entityMutex_[entity];
 129 }
 130
 131
 132 // Return the edge index for a provided edge
 133 inline label dynamicTopoFvMesh::getEdgeIndex
 134 (
 135     const edge& edgeToCheck
 136 ) const
 137 {
 138     if (is2D())
 139     {
 140         // No efficient search method. Loop through all edges.
 141         forAll(edges_, edgeI)
 142         {
 143             if (edges_[edgeI] == edgeToCheck)
 144             {
 145                 return edgeI;
 146             }
 147         }
 148     }
 149     else
 150     {
 151         // Look througg pointEdges list
 152         const labelList& pEdges = pointEdges_[edgeToCheck.start()];
 153
 154         forAll(pEdges, edgeI)
 155         {
 156             if (edges_[pEdges[edgeI]] == edgeToCheck)
 157             {
 158                 return pEdges[edgeI];
 159             }
 160         }
 161     }
 162
 163     // Probably not a valid edge.
 164     FatalErrorIn
 165     (
 166         "inline label dynamicTopoFvMesh::getEdgeIndex"
 167         "(const edge& edgeToCheck) const"
 168     )
 169         << "Could not find an appropriate edge index for edge:"
 170         << edgeToCheck
 171         << abort(FatalError);
 172
 173     return -1;
 174 }
 175
 176
 177 // Given any quad face, pick out a boundary edge that
 178 // contains a triangular face. For 2D simplical meshes only.
 179 inline label dynamicTopoFvMesh::getTriBoundaryEdge
 180 (
 181     const label fIndex
 182 ) const
 183 {
 184     const labelList& fEdges = faceEdges_[fIndex];
 185
 186     forAll(fEdges, edgeI)
 187     {
 188         // Obtain edgeFaces for this edge
 189         const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
 190
 191         forAll(eFaces, faceI)
 192         {
 193             if (faces_[eFaces[faceI]].size() == 3)
 194             {
 195                 // Found a triangular face. Return this edge.
 196                 return fEdges[edgeI];
 197             }
 198         }
 199     }
 200
 201     // This bit should never happen.
 202     FatalErrorIn
 203     (
 204         "inline label dynamicTopoFvMesh::getTriBoundaryEdge"
 205         "(const label fIndex) const"
 206     )
 207         << "Cannot find a triangular face bordering face: "
 208         << fIndex << " :: " << faces_[fIndex]
 209         << abort(FatalError);
 210
 211     return -1;
 212 }
 213
 214
 215 // Check for edge bisection
 216 inline bool dynamicTopoFvMesh::checkBisection
 217 (
 218     const label index
 219 ) const
 220 {
 221     scalar length = 0.0, scale = 0.0;
 222
 223     if (is2D())
 224     {
 225         // If this entity was deleted, skip it.
 226         if (faces_[index].empty())
 227         {
 228             return false;
 229         }
 230
 231         // Fetch the edge
 232         const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
 233
 234         // Measure the boundary edge-length of the face in question
 235         length =
 236         (
 237             linePointRef
 238             (
 239                 points_[edgeToCheck.start()],
 240                 points_[edgeToCheck.end()]
 241             ).mag()
 242         );
 243
 244         // Determine the length-scale at this face
 245         scale = faceLengthScale(index);
 246     }
 247     else
 248     {
 249         // If this entity was deleted, skip it.
 250         if (edgeFaces_[index].empty())
 251         {
 252             return false;
 253         }
 254
 255         // Fetch the edge
 256         const edge& edgeToCheck = edges_[index];
 257
 258         // Measure the edge-length
 259         length =
 260         (
 261             linePointRef
 262             (
 263                 points_[edgeToCheck.start()],
 264                 points_[edgeToCheck.end()]
 265             ).mag()
 266         );
 267
 268         // Determine the length-scale at this point in the mesh
 269         scale = edgeLengthScale(index);
 270     }
 271
 272     if (length > lengthEstimator().ratioMax() * scale)
 273     {
 274         return true;
 275     }
 276     else
 277     {
 278         return false;
 279     }
 280 }
 281
 282
 283 // Check for edge collapse
 284 inline bool dynamicTopoFvMesh::checkCollapse
 285 (
 286     const label index
 287 ) const
 288 {
 289     scalar length = 0.0, scale = 0.0;
 290
 291     if (is2D())
 292     {
 293         // If this entity was deleted, skip it.
 294         if (faces_[index].empty())
 295         {
 296             return false;
 297         }
 298
 299         // Fetch the edge
 300         const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
 301
 302         // Measure the boundary edge-length of the face in question
 303         length =
 304         (
 305             linePointRef
 306             (
 307                 points_[edgeToCheck.start()],
 308                 points_[edgeToCheck.end()]
 309             ).mag()
 310         );
 311
 312         // Determine the length-scale at this face
 313         scale = faceLengthScale(index);
 314     }
 315     else
 316     {
 317         // If this entity was deleted, skip it.
 318         if (edgeFaces_[index].empty())
 319         {
 320             return false;
 321         }
 322
 323         // Fetch the edge
 324         const edge& edgeToCheck = edges_[index];
 325
 326         // Measure the edge-length
 327         length =
 328         (
 329             linePointRef
 330             (
 331                 points_[edgeToCheck.start()],
 332                 points_[edgeToCheck.end()]
 333             ).mag()
 334         );
 335
 336         // Determine the length-scale at this point in the mesh
 337         scale = edgeLengthScale(index);
 338     }
 339
 340     if (length < lengthEstimator().ratioMin() * scale)
 341     {
 342         return true;
 343     }
 344     else
 345     {
 346         return false;
 347     }
 348 }
 349
 350
 351 // Return the entity stack
 352 inline Stack& dynamicTopoFvMesh::stack
 353 (
 354     const label threadID
 355 )
 356 {
 357     return entityStack_[threadID];
 358 }
 359
 360
 361 // Return the integer ID for a given thread
 362 // Return zero for single-threaded operation
 363 inline label dynamicTopoFvMesh::self() const
 364 {
 365     if (threader_->multiThreaded())
 366     {
 367         for (label i = 1; i <= threader_->getNumThreads(); i++)
 368         {
 369             if (handlerPtr_[i].self())
 370             {
 371                 return i;
 372             }
 373         }
 374     }
 375
 376     return 0;
 377 }
 378
 379
 380 // Initialize edge-stacks
 381 inline void dynamicTopoFvMesh::initStacks
 382 (
 383     const labelHashSet& entities
 384 )
 385 {
 386     forAll(entityStack_, stackI)
 387     {
 388         entityStack_[stackI].clear();
 389     }
 390
 391     // Prepare a filling sequence based on threading operation
 392     label tIndex = 0;
 393     labelList tID(threader().getNumThreads());
 394
 395     if (threader_->multiThreaded())
 396     {
 397         forAll(tID, tI)
 398         {
 399             tID[tI] = (tI + 1);
 400         }
 401     }
 402     else
 403     {
 404         tID = 0;
 405     }
 406
 407     if (is2D())
 408     {
 409         forAll(faces_, faceI)
 410         {
 411             // For coupled meshes, avoid certain faces.
 412             if (patchCoupling_.size() || procIndices_.size())
 413             {
 414                 if (entities.found(faceI))
 415                 {
 416                     continue;
 417                 }
 418             }
 419
 420             if (faces_[faceI].size() == 4)
 421             {
 422                 stack(tID[tIndex]).insert(faceI);
 423
 424                 tIndex = tID.fcIndex(tIndex);
 425             }
 426         }
 427     }
 428     else
 429     {
 430         forAll(edges_, edgeI)
 431         {
 432             // For coupled meshes, avoid certain edges.
 433             if (patchCoupling_.size() || procIndices_.size())
 434             {
 435                 if (entities.found(edgeI))
 436                 {
 437                     continue;
 438                 }
 439             }
 440
 441             if (edgeFaces_[edgeI].size())
 442             {
 443                 stack(tID[tIndex]).insert(edgeI);
 444
 445                 tIndex = tID.fcIndex(tIndex);
 446             }
 447         }
 448     }
 449
 450     if (debug > 3 && Pstream::parRun())
 451     {
 452         Pout<< nl << "Stack size: " << stack(0).size() << endl;
 453
 454         if (debug > 4)
 455         {
 456             // Write out stack entities
 457             labelList stackElements(stack(0).size(), -1);
 458
 459             forAll(stackElements, elemI)
 460             {
 461                 stackElements[elemI] = stack(0)[elemI];
 462             }
 463
 464             label elemType = is2D() ? 2 : 1;
 465
 466             writeVTK
 467             (
 468                 "Stack_"
 469               + Foam::name(Pstream::myProcNo()),
 470                 stackElements,
 471                 elemType
 472             );
 473         }
 474     }
 475 }
 476
 477
 478 // Method to determine the old boundary patch index for a given face
 479 // Similar to the polyBoundaryMesh routine, but works on local information
 480 inline label dynamicTopoFvMesh::whichPatch
 481 (
 482     const label index
 483 ) const
 484 {
 485     if (index < nOldInternalFaces_)
 486     {
 487         return -1;
 488     }
 489
 490     for (label i = 0; i < boundaryMesh().size(); i++)
 491     {
 492         if
 493         (
 494             index >= oldPatchStarts_[i]
 495          && index < oldPatchStarts_[i] + oldPatchSizes_[i]
 496         )
 497         {
 498             return i;
 499         }
 500     }
 501
 502     // If not in any of the above, it's possible that the face was added
 503     // at the end of the list. Check addedFacePatches_ for the patch info
 504     Map<label>::const_iterator it = addedFacePatches_.find(index);
 505
 506     if (it != addedFacePatches_.end())
 507     {
 508         return it();
 509     }
 510     else
 511     {
 512         FatalErrorIn
 513         (
 514             "label dynamicTopoFvMesh::whichPatch"
 515             "(const label index) const"
 516         )
 517             << "Cannot find patch information for face index: "
 518             << index << nl
 519             << " It appears that face ordering is"
 520             << " inconsistent with patch information." << nl
 521             << " Mesh info: " << nl
 522             << "    nOldInternalFaces: " << nOldInternalFaces_ << nl
 523             << "    oldPatchStarts: " << oldPatchStarts_ << nl
 524             << "    oldPatchSizes: " << oldPatchSizes_ << nl
 525             << abort(FatalError);
 526     }
 527
 528     return -2;
 529 }
 530
 531
 532 // Method to determine the old boundary patch index for a given edge
 533 inline label dynamicTopoFvMesh::whichEdgePatch
 534 (
 535     const label index
 536 ) const
 537 {
 538     if (index < nOldInternalEdges_)
 539     {
 540         return -1;
 541     }
 542
 543     for (label i = 0; i < boundaryMesh().size(); i++)
 544     {
 545         if
 546         (
 547             index >= oldEdgePatchStarts_[i]
 548          && index < oldEdgePatchStarts_[i] + oldEdgePatchSizes_[i]
 549         )
 550         {
 551             return i;
 552         }
 553     }
 554
 555     // If not in any of the above, it's possible that the edge was added
 556     // at the end of the list. Check addedEdgePatches_ for the patch info
 557     Map<label>::const_iterator it = addedEdgePatches_.find(index);
 558
 559     if (it != addedEdgePatches_.end())
 560     {
 561         return it();
 562     }
 563     else
 564     {
 565         FatalErrorIn
 566         (
 567             "label dynamicTopoFvMesh::whichEdgePatch"
 568             "(const label index) const"
 569         )
 570             << "Cannot find patch information for edge index: "
 571             << index << nl
 572             << " It appears that edge ordering is"
 573             << " inconsistent with patch information." << nl
 574             << " Mesh info: " << nl
 575             << "    nOldInternalEdges: " << nOldInternalEdges_ << nl
 576             << "    oldEdgePatchStarts: " << oldEdgePatchStarts_ << nl
 577             << "    oldEdgePatchSizes: " << oldEdgePatchSizes_ << nl
 578             << abort(FatalError);
 579     }
 580
 581     return -2;
 582 }
 583
 584
 585 // Report or alter topo-modification status
 586 inline label& dynamicTopoFvMesh::status(const label type)
 587 {
 588     if (type < 0 || type >= TOTAL_OP_TYPES)
 589     {
 590         FatalErrorIn
 591         (
 592             "inline label& dynamicTopoFvMesh::status(const label type)"
 593         )
 594             << " Unknown type: " << type << nl
 595             << abort(FatalError);
 596     }
 597
 598     return statistics_[type];
 599 }
 600
 601
 602 // Set a particular face index as flipped.
 603 inline void dynamicTopoFvMesh::setFlip(const label fIndex)
 604 {
 605     if (fIndex < nOldFaces_)
 606     {
 607         labelHashSet::iterator it = flipFaces_.find(fIndex);
 608
 609         if (it == flipFaces_.end())
 610         {
 611             flipFaces_.insert(fIndex);
 612         }
 613         else
 614         {
 615             flipFaces_.erase(it);
 616         }
 617     }
 618 }
 619
 620
 621 // Check for processor priority
 622 template <class BinaryOp>
 623 inline bool dynamicTopoFvMesh::priority
 624 (
 625     const label first,
 626     const BinaryOp& bop,
 627     const label second
 628 ) const
 629 {
 630     return bop(procPriority_[first], procPriority_[second]);
 631 }
 632
 633
 634 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 635
 636 } // End namespace Foam
 637
 638 // ************************************************************************* //