src/mesh/cfMesh/utilities/helperFunctions/helperFunctionsTopologyManipulationI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "error.H"
  29 #include "helperFunctionsTopologyManipulation.H"
  30 #include "edgeList.H"
  31 #include "pointField.H"
  32 #include "boolList.H"
  33 #include "cellList.H"
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  39
  40 namespace help
  41 {
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  44
  45 template<class faceType1, class faceType2>
  46 bool areFacesEqual(const faceType1& f1, const faceType2& f2)
  47 {
  48     //- their size mut be equal
  49     if( f1.size() != f2.size() )
  50         return false;
  51
  52     //- find the starting point the the second face
  53     label start(-1);
  54     const label s = f2.size();
  55     bool equalOrientation(false);
  56
  57     forAll(f2, pI)
  58     {
  59         if( f1[0] == f2[pI] )
  60         {
  61             start = pI;
  62
  63             if( f1[1] == f2[(pI+1)%s] )
  64             {
  65                 //- faces have the same orientation
  66                 equalOrientation = true;
  67             }
  68             else if( f1[1] != f2[(s-1+pI)%s] )
  69             {
  70                 //- faces shall have opposite orientation, but do not match
  71                 return false;
  72             }
  73         }
  74     }
  75
  76     if( start < 0 )
  77         return false;
  78
  79     if( equalOrientation )
  80     {
  81         //- same orientation, check if all points match
  82         for(label pI=1;pI<s;++pI)
  83         {
  84             if( f1[pI] != f2[(pI+start)%s] )
  85                 return false;
  86         }
  87     }
  88     else
  89     {
  90         //- faces with opposite orientation, check if all points match
  91         for(label pI=1;pI<s;++pI)
  92         {
  93             if( f1[pI] != f2[(start+s-pI)%s] )
  94                 return false;
  95         }
  96     }
  97
  98     //- faces are equal
  99     return true;
 100 }
 101
 102 template<class T, class ListType >
 103 label positionInList(const T& elmt, const ListType& l)
 104 {
 105     for(label i=0;i<l.size();++i)
 106         if( l[i] == elmt )
 107             return i;
 108
 109     return -1;
 110 }
 111
 112 template<class faceType>
 113 faceType reverseFace(const faceType& f)
 114 {
 115     faceType rf;
 116     rf.setSize(f.size());
 117
 118     rf[0] = f[0];
 119
 120     const label size = f.size();
 121     for(label i=1;i<size;++i)
 122         rf[f.size()-i] = f[i];
 123
 124     return rf;
 125 }
 126
 127 template<class faceType1, class faceType2>
 128 inline face mergeTwoFaces(const faceType1& f1, const faceType2& f2)
 129 {
 130     DynList<bool> ce1, ce2;
 131     ce1.setSize(f1.size());
 132     ce1 = false;
 133     ce2.setSize(f2.size());
 134     ce2 = false;
 135
 136     forAll(f1, eI)
 137     {
 138         const edge e1(f1[eI], f1[f1.fcIndex(eI)]);
 139
 140         forAll(f2, eJ)
 141         {
 142             const edge e2(f2[eJ], f2[f2.fcIndex(eJ)]);
 143
 144             if( e1 == e2 )
 145             {
 146                 ce1[eI] = true;
 147                 ce2[eJ] = true;
 148                 break;
 149             }
 150         }
 151     }
 152
 153     DynList<edge> fEdges;
 154     forAll(ce1, eI)
 155     {
 156         if( !ce1[eI] )
 157         {
 158             const edge e1(f1[eI], f1[f1.fcIndex(eI)]);
 159             fEdges.append(e1);
 160         }
 161     }
 162
 163     forAll(ce2, eI)
 164     {
 165         if( !ce2[eI] )
 166         {
 167             const edge e2(f2[eI], f2[f2.fcIndex(eI)]);
 168             fEdges.append(e2);
 169         }
 170     }
 171
 172     return face(sortEdgeChain(fEdges));
 173 }
 174
 175 inline edgeList modifyFacesToShareOneEdge(face& f1, face& f2)
 176 {
 177     const edgeList edg1 = f1.edges();
 178     const edgeList edg2 = f2.edges();
 179
 180     boolList she1(f1.size(), false);
 181     boolList she2(f2.size(), false);
 182
 183     edgeList sharedEdges(edg1);
 184     label nSharedEdges(0);
 185
 186     forAll(edg1, eI)
 187         forAll(edg2, eJ)
 188             if( edg1[eI] == edg2[eJ] )
 189             {
 190                 sharedEdges.newElmt(nSharedEdges++) = edg1[eI];
 191
 192                 she1[eI] = true;
 193                 she2[eJ] = true;
 194                 break;
 195             }
 196
 197     face newF(f1);
 198     label i(0);
 199     forAll(f1, pI)
 200         if( !(she1[pI] && she1[(pI-1+f1.size())%f1.size()]) )
 201             newF[i++] = f1[pI];
 202
 203     newF.setSize(i);
 204     f1 = newF;
 205
 206     newF = f2;
 207     i = 0;
 208     forAll(f2, pI)
 209         if( !(she2[pI] && she2[(pI-1+f2.size())%f2.size()]) )
 210             newF[i++] = f2[pI];
 211
 212     newF.setSize(i);
 213     f2 = newF;
 214
 215     sharedEdges.setSize(nSharedEdges);
 216     return sharedEdges;
 217 }
 218
 219 inline face createFaceFromRemovedPart(const face& fOrig, const face& fCut)
 220 {
 221     if( fCut.size() == 0 )
 222         return fOrig;
 223
 224     const edgeList eOrig = fOrig.edges();
 225     const edgeList eCut = fCut.edges();
 226
 227     boolList usedEdge(eOrig.size(), false);
 228
 229     forAll(eOrig, eI)
 230         forAll(eCut, eJ)
 231             if( eOrig[eI] == eCut[eJ] )
 232             {
 233                 usedEdge[eI] = true;
 234                 break;
 235             }
 236
 237     face f(fOrig);
 238     direction i(0);
 239
 240     forAll(fOrig, pI)
 241         if( !(usedEdge[pI] && usedEdge[(pI-1+fOrig.size())%fOrig.size()]) )
 242         {
 243             f[i++] = fOrig[pI];
 244         }
 245
 246     f.setSize(i);
 247
 248     return f;
 249 }
 250
 251 inline face removeEdgesFromFace
 252 (
 253     const face& fOrig,
 254     const DynList<edge>& removeEdges
 255 )
 256 {
 257     boolList foundEdge(fOrig.size(), false);
 258
 259     forAll(removeEdges, reI)
 260         forAll(fOrig, eI)
 261             if( removeEdges[reI] == fOrig.faceEdge(eI) )
 262             {
 263                 foundEdge[eI] = true;
 264                 break;
 265             }
 266
 267     face newF(fOrig.size());
 268     label i(0);
 269
 270     forAll(fOrig, pI)
 271         if( !(foundEdge[pI] && foundEdge[fOrig.rcIndex(pI)]) )
 272             newF[i++] = fOrig[pI];
 273
 274     newF.setSize(i);
 275
 276     return newF;
 277 }
 278
 279 inline void findOpenEdges(const faceList& cellFaces, DynList<edge>& openEdges)
 280 {
 281     DynList<edge> cellEdges;
 282     DynList<label> nAppearances;
 283
 284     forAll(cellFaces, fI)
 285     {
 286         const edgeList edges = cellFaces[fI].edges();
 287
 288         forAll(edges, eI)
 289         {
 290             const label pos = cellEdges.containsAtPosition(edges[eI]);
 291
 292             if( pos == -1 )
 293             {
 294                 cellEdges.append(edges[eI]);
 295                 nAppearances.append(1);
 296             }
 297             else
 298             {
 299                 nAppearances[pos]++;
 300             }
 301         }
 302     }
 303
 304     openEdges.setSize(12);
 305     openEdges.clear();
 306
 307     forAll(nAppearances, eI)
 308         if( nAppearances[eI] == 1 )
 309         {
 310             openEdges.append(cellEdges[eI]);
 311         }
 312         else if( nAppearances[eI] > 2 )
 313         {
 314             FatalErrorIn
 315             (
 316                 "void findOpenEdges(const faceList& cellFaces,"
 317                 "DynList<edge>& openEdges)"
 318             ) << "More than two faces in " << cellFaces
 319                 << " share edge " << cellEdges[eI] << abort(FatalError);
 320         }
 321 }
 322
 323 template<class faceType1, class faceType2>
 324 inline bool shareAnEdge(const faceType1& f1, const faceType2& f2)
 325 {
 326     forAll(f1, eI)
 327     {
 328         const edge e1(f1[eI], f1[f1.fcIndex(eI)]);
 329
 330         forAll(f2, eJ)
 331         {
 332             const edge e2(f2[eJ], f2[f2.fcIndex(eJ)]);
 333
 334             if( e1 == e2 )
 335                 return true;
 336         }
 337     }
 338
 339     return false;
 340 }
 341
 342 template<class faceType1, class faceType2>
 343 inline edge sharedEdge(const faceType1& f1, const faceType2& f2)
 344 {
 345     forAll(f1, eI)
 346     {
 347         const edge e1(f1[eI], f1[f1.fcIndex(eI)]);
 348
 349         forAll(f2, eJ)
 350         {
 351             const edge e2(f2[eJ], f2[f2.fcIndex(eJ)]);
 352
 353             if( e1 == e2 )
 354                 return e1;
 355         }
 356     }
 357
 358     return edge(-1, -1);
 359 }
 360
 361 template<class faceType>
 362 inline label positionOfEdgeInFace(const edge& e, const faceType& f)
 363 {
 364     forAll(f, eI)
 365     {
 366         const edge fe(f[eI], f[f.fcIndex(eI)]);
 367
 368         if( fe == e )
 369             return eI;
 370     }
 371
 372     return -1;
 373 }
 374
 375 template<class faceType1, class faceType2>
 376 inline label sharedVertex(const faceType1& f1, const faceType2& f2)
 377 {
 378     forAll(f1, pI)
 379         forAll(f2, pJ)
 380             if( f1[pI] == f2[pJ] )
 381                 return f1[pI];
 382
 383     return -1;
 384 }
 385
 386 template<class faceType1, class faceType2>
 387 inline bool shareAVertex(const faceType1& f1, const faceType2& f2)
 388 {
 389     forAll(f1, pI)
 390         forAll(f2, pJ)
 391             if( f1[pI] == f2[pJ] )
 392                 return true;
 393
 394     return false;
 395 }
 396
 397 template<class faceListType>
 398 inline label sharedVertex(const faceListType& fcs)
 399 {
 400     forAll(fcs[0], pI)
 401     {
 402         bool allFound(true);
 403
 404         for(label i=1;i<fcs.size();++i)
 405         {
 406             bool found(false);
 407             forAll(fcs[i], pJ)
 408                 if( fcs[0][pI] == fcs[i][pJ] )
 409                     found = true;
 410
 411             if( !found )
 412             {
 413                 allFound = false;
 414                 break;
 415             }
 416         }
 417
 418         if( allFound )
 419             return fcs[0][pI];
 420     }
 421
 422     return -1;
 423 }
 424
 425 template<class boolListType>
 426 inline bool areElementsInChain(const boolListType& sel)
 427 {
 428     DynList<bool> selInChain(sel.size(), false);
 429
 430     forAll(sel, eI)
 431     {
 432         if( sel[eI] )
 433         {
 434             selInChain[eI] = true;
 435             bool found;
 436             do
 437             {
 438                 found = false;
 439                 forAll(selInChain, eJ)
 440                     if(
 441                         !selInChain[eJ] && sel[eJ] &&
 442                         (
 443                             selInChain[sel.fcIndex(eJ)] ||
 444                             selInChain[sel.rcIndex(eJ)]
 445                         )
 446                     )
 447                     {
 448                         found = true;
 449                         selInChain[eJ] = true;
 450                     }
 451             } while( found );
 452
 453             break;
 454         }
 455     }
 456
 457     forAll(sel, eI)
 458     {
 459         if( sel[eI] && !selInChain[eI] )
 460             return false;
 461     }
 462
 463     return true;
 464 }
 465
 466 inline void zipOpenChain(DynList<edge>& bEdges)
 467 {
 468     //- close the chain if open
 469     DynList<label> chainVertices;
 470     List<label> nAppearances;
 471     forAll(bEdges, eI)
 472     {
 473         const edge& e = bEdges[eI];
 474         forAll(e, pI)
 475         {
 476             const label pos = chainVertices.containsAtPosition(e[pI]);
 477
 478             if( pos == -1 )
 479             {
 480                 nAppearances[chainVertices.size()] = 1;
 481                 chainVertices.append(e[pI]);
 482             }
 483             else
 484             {
 485                 ++nAppearances[pos];
 486             }
 487         }
 488     }
 489
 490     bool closed(true);
 491     DynList<label> openVertices(2);
 492     forAll(chainVertices, pI)
 493         if( nAppearances[pI] == 1 )
 494         {
 495             closed = false;
 496             openVertices.append(chainVertices[pI]);
 497         }
 498
 499     if( !closed && (openVertices.size() == 2) )
 500     {
 501         forAll(bEdges, eI)
 502             if( bEdges[eI].end() == openVertices[0] )
 503             {
 504                 bEdges.append(edge(openVertices[0], openVertices[1]));
 505                 break;
 506             }
 507             else if( bEdges[eI].end() == openVertices[1] )
 508             {
 509                 bEdges.append(edge(openVertices[1], openVertices[0]));
 510                 break;
 511             }
 512     }
 513     else if( !closed )
 514     {
 515         FatalErrorIn
 516         (
 517             "void dualMeshExtractor::decomposeCreatedPoly::"
 518             "createMissingFaces(List<faceList>& cFaces)"
 519         ) << "Chain has " << openVertices << " open vertices"
 520             << abort(FatalError);
 521     }
 522 }
 523
 524 inline labelList sortEdgeChain(const DynList<edge>& bEdges)
 525 {
 526     boolList sorted(bEdges.size(), false);
 527
 528     DynList<edge> sortedEdges;
 529     sortedEdges.append(bEdges[0]);
 530     sorted[0] = true;
 531     direction i(0);
 532
 533     bool finished;
 534     do
 535     {
 536         finished = true;
 537
 538         forAll(bEdges, eI)
 539             if( !sorted[eI] )
 540             {
 541                 if( sortedEdges[i].end() == bEdges[eI].start() )
 542                 {
 543                     sorted[eI] = true;
 544                     finished = false;
 545                     sortedEdges.append(bEdges[eI]);
 546                     ++i;
 547                 }
 548                 else if( sortedEdges[i].end() == bEdges[eI].end() )
 549                 {
 550                     FatalErrorIn
 551                     (
 552                         "labelList sortEdgeChain("
 553                         "const DynList<edge>& bEdges)"
 554                     ) << "Chain is not oriented correctly!"
 555                         << abort(FatalError);
 556                 }
 557             }
 558     } while( !finished );
 559
 560     labelList sortPoints(bEdges.size());
 561     forAll(sortedEdges, eI)
 562         sortPoints[eI] = sortedEdges[eI].start();
 563
 564     return sortPoints;
 565 }
 566
 567 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 568
 569 } // End namespace help
 570
 571 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 572
 573 } // End namespace Foam
 574
 575 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //