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-2008 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 "labelHashSet.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 = (fp + 1) % f.size();
  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 = (face0Start + 1) % f0.size();
 150
 151     const face& f1 = mesh_.faces()[face1];
 152
 153     label face1Start = findIndex(f1, e.start());
 154     label face1End = (face1Start + 1) % f1.size();
 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].size() == 0)
 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                     superFace.shrink();
 373
 374                     faces[superFaceI] = face(superFace);
 375                 }
 376             }
 377         }
 378     }
 379 }
 380
 381
 382 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 383
 384 // Construct from components
 385 Foam::cellFeatures::cellFeatures
 386 (
 387     const primitiveMesh& mesh,
 388     const scalar minCos,
 389     const label cellI
 390 )
 391 :
 392     mesh_(mesh),
 393     minCos_(minCos),
 394     cellI_(cellI),
 395     featureEdge_(10*mesh.cellEdges()[cellI].size()),
 396     facesPtr_(NULL),
 397     faceMap_(0)
 398 {
 399     const labelList& cEdges = mesh_.cellEdges()[cellI_];
 400
 401     forAll(cEdges, cEdgeI)
 402     {
 403         label edgeI = cEdges[cEdgeI];
 404
 405         if (isCellFeatureEdge(minCos_, edgeI))
 406         {
 407             featureEdge_.insert(edgeI);
 408         }
 409     }
 410
 411 }
 412
 413
 414 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 415
 416 Foam::cellFeatures::~cellFeatures()
 417 {
 418     deleteDemandDrivenData(facesPtr_);
 419 }
 420
 421
 422 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 423
 424 bool Foam::cellFeatures::isFeaturePoint(const label edge0, const label edge1)
 425  const
 426 {
 427     if
 428     (
 429         (edge0 < 0)
 430      || (edge0 >= mesh_.nEdges())
 431      || (edge1 < 0)
 432      || (edge1 >= mesh_.nEdges())
 433     )
 434     {
 435         FatalErrorIn
 436         (
 437             "cellFeatures::isFeatureVertex(const label, const label)"
 438         )   << "Illegal edge labels : edge0:" << edge0 << " edge1:" << edge1
 439             << abort(FatalError);
 440     }
 441
 442     const edge& e0 = mesh_.edges()[edge0];
 443
 444     vector e0Vec = e0.vec(mesh_.points());
 445     e0Vec /= mag(e0Vec);
 446
 447     const edge& e1 = mesh_.edges()[edge1];
 448
 449     vector e1Vec = e1.vec(mesh_.points());
 450     e1Vec /= mag(e1Vec);
 451
 452     scalar cosAngle;
 453
 454     if
 455     (
 456         (e0.start() == e1.end())
 457      || (e0.end() == e1.start())
 458     )
 459     {
 460         // Same direction
 461         cosAngle = e0Vec & e1Vec;
 462     }
 463     else if
 464     (
 465         (e0.start() == e1.start())
 466      || (e0.end() == e1.end())
 467     )
 468     {
 469         // back on back
 470         cosAngle = - e0Vec & e1Vec;
 471     }
 472     else
 473     {
 474         cosAngle = GREAT;   // satisfy compiler
 475
 476         FatalErrorIn
 477         (
 478             "cellFeatures::isFeaturePoint(const label, const label"
 479             ", const label)"
 480         )   << "Edges do not share common vertex. e0:" << e0
 481             << " e1:" << e1 << abort(FatalError);
 482     }
 483
 484     if (cosAngle < minCos_)
 485     {
 486         // Angle larger than criterium
 487         return true;
 488     }
 489     else
 490     {
 491         return false;
 492     }
 493 }
 494
 495
 496 bool Foam::cellFeatures::isFeatureVertex(const label faceI, const label vertI)
 497  const
 498 {
 499     if
 500     (
 501         (faceI < 0)
 502      || (faceI >= mesh_.nFaces())
 503      || (vertI < 0)
 504      || (vertI >= mesh_.nPoints())
 505     )
 506     {
 507         FatalErrorIn
 508         (
 509             "cellFeatures::isFeatureVertex(const label, const label)"
 510         )   << "Illegal face " << faceI << " or vertex " << vertI
 511             << abort(FatalError);
 512     }
 513
 514     const labelList& pEdges = mesh_.pointEdges()[vertI];
 515
 516     label edge0 = -1;
 517     label edge1 = -1;
 518
 519     forAll(pEdges, pEdgeI)
 520     {
 521         label edgeI = pEdges[pEdgeI];
 522
 523         if (meshTools::edgeOnFace(mesh_, faceI, edgeI))
 524         {
 525             if (edge0 == -1)
 526             {
 527                 edge0 = edgeI;
 528             }
 529             else
 530             {
 531                 edge1 = edgeI;
 532
 533                 // Found the two edges.
 534                 break;
 535             }
 536         }
 537     }
 538
 539     if (edge1 == -1)
 540     {
 541         FatalErrorIn
 542         (
 543             "cellFeatures::isFeatureVertex(const label, const label)"
 544         )   << "Did not find two edges sharing vertex " << vertI
 545             << " on face " << faceI << " vertices:" << mesh_.faces()[faceI]
 546             << abort(FatalError);
 547     }
 548
 549     return isFeaturePoint(edge0, edge1);
 550 }
 551
 552
 553 // ************************************************************************* //