src/meshTools/cellFeatures/cellFeatures.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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "cellFeatures.H"
  29 #include "primitiveMesh.H"
  30 #include "HashSet.H"
  31 #include "Map.H"
  32 #include "demandDrivenData.H"
  33 #include "ListOps.H"
  34 #include "meshTools.H"
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38
  39 // Return true if edge start and end are on increasing face vertices. (edge is
  40 // guaranteed to be on face)
  41 bool Foam::cellFeatures::faceAlignedEdge(const label faceI, const label edgeI)
  42  const
  43 {
  44     const edge& e = mesh_.edges()[edgeI];
  45
  46     const face& f = mesh_.faces()[faceI];
  47
  48     forAll(f, fp)
  49     {
  50         if (f[fp] == e.start())
  51         {
  52             label fp1 = f.fcIndex(fp);
  53
  54             return f[fp1] == e.end();
  55         }
  56     }
  57
  58     FatalErrorIn
  59     (
  60         "cellFeatures::faceAlignedEdge(const label, const label)"
  61     )   << "Can not find edge " << mesh_.edges()[edgeI]
  62         << " on face " << faceI << abort(FatalError);
  63
  64     return false;
  65 }
  66
  67
  68 // Return edge in featureEdge that uses vertI and is on same superface
  69 // but is not edgeI
  70 Foam::label Foam::cellFeatures::nextEdge
  71 (
  72     const Map<label>& toSuperFace,
  73     const label superFaceI,
  74     const label thisEdgeI,
  75     const label thisVertI
  76 ) const
  77 {
  78     const labelList& pEdges = mesh_.pointEdges()[thisVertI];
  79
  80     forAll(pEdges, pEdgeI)
  81     {
  82         label edgeI = pEdges[pEdgeI];
  83
  84         if ((edgeI != thisEdgeI) && featureEdge_.found(edgeI))
  85         {
  86             // Check that edge is used by a face on same superFace
  87
  88             const labelList& eFaces = mesh_.edgeFaces()[edgeI];
  89
  90             forAll(eFaces, eFaceI)
  91             {
  92                 label faceI = eFaces[eFaceI];
  93
  94                 if
  95                 (
  96                     meshTools::faceOnCell(mesh_, cellI_, faceI)
  97                  && (toSuperFace[faceI] == superFaceI)
  98                 )
  99                 {
 100                     return edgeI;
 101                 }
 102             }
 103         }
 104     }
 105
 106     FatalErrorIn
 107     (
 108         "cellFeatures::nextEdge(const label, const Map<label>"
 109         ", const labelHashSet&, const label, const label, const label)"
 110     )   << "Can not find edge in " << featureEdge_ << " connected to edge "
 111         << thisEdgeI << " at vertex " << thisVertI << endl
 112         << "This might mean that the externalEdges do not form a closed loop"
 113         << abort(FatalError);
 114
 115     return -1;
 116 }
 117
 118
 119 // Return true if angle between faces using it is larger than certain value.
 120 bool Foam::cellFeatures::isCellFeatureEdge
 121 (
 122     const scalar minCos,
 123     const label edgeI
 124 ) const
 125 {
 126     // Get the two faces using this edge
 127
 128     label face0;
 129     label face1;
 130     meshTools::getEdgeFaces(mesh_, cellI_, edgeI, face0, face1);
 131
 132     // Check the angle between them by comparing the face normals.
 133
 134     vector n0 = mesh_.faceAreas()[face0];
 135     n0 /= mag(n0);
 136
 137     vector n1 = mesh_.faceAreas()[face1];
 138     n1 /= mag(n1);
 139
 140     scalar cosAngle = n0 & n1;
 141
 142
 143     const edge& e = mesh_.edges()[edgeI];
 144
 145     const face& f0 = mesh_.faces()[face0];
 146
 147     label face0Start = findIndex(f0, e.start());
 148     label face0End   = f0.fcIndex(face0Start);
 149
 150     const face& f1 = mesh_.faces()[face1];
 151
 152     label face1Start = findIndex(f1, e.start());
 153     label face1End   = f1.fcIndex(face1Start);
 154
 155     if
 156     (
 157         (
 158             (f0[face0End] == e.end())
 159          && (f1[face1End] != e.end())
 160         )
 161      || (
 162             (f0[face0End] != e.end())
 163          && (f1[face1End] == e.end())
 164         )
 165     )
 166     {
 167     }
 168     else
 169     {
 170         cosAngle = -cosAngle;
 171     }
 172
 173     if (cosAngle < minCos)
 174     {
 175         return true;
 176     }
 177     else
 178     {
 179         return false;
 180     }
 181 }
 182
 183
 184 // Recursively mark (on toSuperFace) all face reachable across non-feature
 185 // edges.
 186 void Foam::cellFeatures::walkSuperFace
 187 (
 188     const label faceI,
 189     const label superFaceI,
 190     Map<label>& toSuperFace
 191 ) const
 192 {
 193     if (!toSuperFace.found(faceI))
 194     {
 195         toSuperFace.insert(faceI, superFaceI);
 196
 197         const labelList& fEdges = mesh_.faceEdges()[faceI];
 198
 199         forAll(fEdges, fEdgeI)
 200         {
 201             label edgeI = fEdges[fEdgeI];
 202
 203             if (!featureEdge_.found(edgeI))
 204             {
 205                 label face0;
 206                 label face1;
 207                 meshTools::getEdgeFaces(mesh_, cellI_, edgeI, face0, face1);
 208
 209                 if (face0 == faceI)
 210                 {
 211                     face0 = face1;
 212                 }
 213
 214                 walkSuperFace
 215                 (
 216                     face0,
 217                     superFaceI,
 218                     toSuperFace
 219                 );
 220             }
 221         }
 222     }
 223 }
 224
 225
 226 void Foam::cellFeatures::calcSuperFaces() const
 227 {
 228     // Determine superfaces by edge walking across non-feature edges
 229
 230     const labelList& cFaces = mesh_.cells()[cellI_];
 231
 232     // Mapping from old to super face:
 233     //    <not found> : not visited
 234     //    >=0 : superFace
 235     Map<label> toSuperFace(10*cFaces.size());
 236
 237     label superFaceI = 0;
 238
 239     forAll(cFaces, cFaceI)
 240     {
 241         label faceI = cFaces[cFaceI];
 242
 243         if (!toSuperFace.found(faceI))
 244         {
 245             walkSuperFace
 246             (
 247                 faceI,
 248                 superFaceI,
 249                 toSuperFace
 250             );
 251             superFaceI++;
 252         }
 253     }
 254
 255     // Construct superFace-to-oldface mapping.
 256
 257     faceMap_.setSize(superFaceI);
 258
 259     forAll(cFaces, cFaceI)
 260     {
 261         label faceI = cFaces[cFaceI];
 262
 263         faceMap_[toSuperFace[faceI]].append(faceI);
 264     }
 265
 266     forAll(faceMap_, superI)
 267     {
 268         faceMap_[superI].shrink();
 269     }
 270
 271
 272     // Construct superFaces
 273
 274     facesPtr_ = new faceList(superFaceI);
 275
 276     faceList& faces = *facesPtr_;
 277
 278     forAll(cFaces, cFaceI)
 279     {
 280         label faceI = cFaces[cFaceI];
 281
 282         label superFaceI = toSuperFace[faceI];
 283
 284         if (faces[superFaceI].empty())
 285         {
 286             // Superface not yet constructed.
 287
 288             // Find starting feature edge on face.
 289             label startEdgeI = -1;
 290
 291             const labelList& fEdges = mesh_.faceEdges()[faceI];
 292
 293             forAll(fEdges, fEdgeI)
 294             {
 295                 label edgeI = fEdges[fEdgeI];
 296
 297                 if (featureEdge_.found(edgeI))
 298                 {
 299                     startEdgeI = edgeI;
 300
 301                     break;
 302                 }
 303             }
 304
 305
 306             if (startEdgeI != -1)
 307             {
 308                 // Walk point-edge-point along feature edges
 309
 310                 dynamicLabelList superFace(10*mesh_.faces()[faceI].size());
 311
 312                 const edge& e = mesh_.edges()[startEdgeI];
 313
 314                 // Walk either start-end or end-start depending on orientation
 315                 // of face. SuperFace will have cellI as owner.
 316                 bool flipOrientation =
 317                     (mesh_.faceOwner()[faceI] == cellI_)
 318                   ^ (faceAlignedEdge(faceI, startEdgeI));
 319
 320                 label startVertI = -1;
 321
 322                 if (flipOrientation)
 323                 {
 324                     startVertI = e.end();
 325                 }
 326                 else
 327                 {
 328                     startVertI = e.start();
 329                 }
 330
 331                 label edgeI = startEdgeI;
 332
 333                 label vertI = e.otherVertex(startVertI);
 334
 335                 do
 336                 {
 337                     label newEdgeI = nextEdge
 338                     (
 339                         toSuperFace,
 340                         superFaceI,
 341                         edgeI,
 342                         vertI
 343                     );
 344
 345                     // Determine angle between edges.
 346                     if (isFeaturePoint(edgeI, newEdgeI))
 347                     {
 348                         superFace.append(vertI);
 349                     }
 350
 351                     edgeI = newEdgeI;
 352
 353                     if (vertI == startVertI)
 354                     {
 355                         break;
 356                     }
 357
 358                     vertI = mesh_.edges()[edgeI].otherVertex(vertI);
 359                 }
 360                 while (true);
 361
 362                 if (superFace.size() <= 2)
 363                 {
 364                     WarningIn("cellFeatures::calcSuperFaces")
 365                         << " Can not collapse faces " << faceMap_[superFaceI]
 366                         << " into one big face on cell " << cellI_ << endl
 367                         << "Try decreasing minCos:" << minCos_ << endl;
 368                 }
 369                 else
 370                 {
 371                     faces[superFaceI].transfer(superFace);
 372                 }
 373             }
 374         }
 375     }
 376 }
 377
 378
 379 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 380
 381 // Construct from components
 382 Foam::cellFeatures::cellFeatures
 383 (
 384     const primitiveMesh& mesh,
 385     const scalar minCos,
 386     const label cellI
 387 )
 388 :
 389     mesh_(mesh),
 390     minCos_(minCos),
 391     cellI_(cellI),
 392     featureEdge_(10*mesh.cellEdges()[cellI].size()),
 393     facesPtr_(NULL),
 394     faceMap_(0)
 395 {
 396     const labelList& cEdges = mesh_.cellEdges()[cellI_];
 397
 398     forAll(cEdges, cEdgeI)
 399     {
 400         label edgeI = cEdges[cEdgeI];
 401
 402         if (isCellFeatureEdge(minCos_, edgeI))
 403         {
 404             featureEdge_.insert(edgeI);
 405         }
 406     }
 407
 408 }
 409
 410
 411 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 412
 413 Foam::cellFeatures::~cellFeatures()
 414 {
 415     deleteDemandDrivenData(facesPtr_);
 416 }
 417
 418
 419 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 420
 421 bool Foam::cellFeatures::isFeaturePoint(const label edge0, const label edge1)
 422  const
 423 {
 424     if
 425     (
 426         (edge0 < 0)
 427      || (edge0 >= mesh_.nEdges())
 428      || (edge1 < 0)
 429      || (edge1 >= mesh_.nEdges())
 430     )
 431     {
 432         FatalErrorIn
 433         (
 434             "cellFeatures::isFeatureVertex(const label, const label)"
 435         )   << "Illegal edge labels : edge0:" << edge0 << " edge1:" << edge1
 436             << abort(FatalError);
 437     }
 438
 439     const edge& e0 = mesh_.edges()[edge0];
 440
 441     vector e0Vec = e0.vec(mesh_.points());
 442     e0Vec /= mag(e0Vec);
 443
 444     const edge& e1 = mesh_.edges()[edge1];
 445
 446     vector e1Vec = e1.vec(mesh_.points());
 447     e1Vec /= mag(e1Vec);
 448
 449     scalar cosAngle;
 450
 451     if
 452     (
 453         (e0.start() == e1.end())
 454      || (e0.end() == e1.start())
 455     )
 456     {
 457         // Same direction
 458         cosAngle = e0Vec & e1Vec;
 459     }
 460     else if
 461     (
 462         (e0.start() == e1.start())
 463      || (e0.end() == e1.end())
 464     )
 465     {
 466         // back on back
 467         cosAngle = - e0Vec & e1Vec;
 468     }
 469     else
 470     {
 471         cosAngle = GREAT;   // satisfy compiler
 472
 473         FatalErrorIn
 474         (
 475             "cellFeatures::isFeaturePoint(const label, const label"
 476             ", const label)"
 477         )   << "Edges do not share common vertex. e0:" << e0
 478             << " e1:" << e1 << abort(FatalError);
 479     }
 480
 481     if (cosAngle < minCos_)
 482     {
 483         // Angle larger than criterium
 484         return true;
 485     }
 486     else
 487     {
 488         return false;
 489     }
 490 }
 491
 492
 493 bool Foam::cellFeatures::isFeatureVertex(const label faceI, const label vertI)
 494  const
 495 {
 496     if
 497     (
 498         (faceI < 0)
 499      || (faceI >= mesh_.nFaces())
 500      || (vertI < 0)
 501      || (vertI >= mesh_.nPoints())
 502     )
 503     {
 504         FatalErrorIn
 505         (
 506             "cellFeatures::isFeatureVertex(const label, const label)"
 507         )   << "Illegal face " << faceI << " or vertex " << vertI
 508             << abort(FatalError);
 509     }
 510
 511     const labelList& pEdges = mesh_.pointEdges()[vertI];
 512
 513     label edge0 = -1;
 514     label edge1 = -1;
 515
 516     forAll(pEdges, pEdgeI)
 517     {
 518         label edgeI = pEdges[pEdgeI];
 519
 520         if (meshTools::edgeOnFace(mesh_, faceI, edgeI))
 521         {
 522             if (edge0 == -1)
 523             {
 524                 edge0 = edgeI;
 525             }
 526             else
 527             {
 528                 edge1 = edgeI;
 529
 530                 // Found the two edges.
 531                 break;
 532             }
 533         }
 534     }
 535
 536     if (edge1 == -1)
 537     {
 538         FatalErrorIn
 539         (
 540             "cellFeatures::isFeatureVertex(const label, const label)"
 541         )   << "Did not find two edges sharing vertex " << vertI
 542             << " on face " << faceI << " vertices:" << mesh_.faces()[faceI]
 543             << abort(FatalError);
 544     }
 545
 546     return isFeaturePoint(edge0, edge1);
 547 }
 548
 549
 550 // ************************************************************************* //