src/meshTools/cellFeatures/cellFeatures.C

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