src/meshTools/cellFeatures/cellFeatures.C

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