src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "primitiveMesh.H"
  28 #include "DynamicList.H"
  29 #include "demandDrivenData.H"
  30 #include "SortableList.H"
  31 #include "ListOps.H"
  32
  33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  34
  35 // Returns edgeI between two points.
  36 Foam::label Foam::primitiveMesh::getEdge
  37 (
  38     List<DynamicList<label> >& pe,
  39     DynamicList<edge>& es,
  40
  41     const label pointI,
  42     const label nextPointI
  43 )
  44 {
  45     // Find connection between pointI and nextPointI
  46     forAll(pe[pointI], ppI)
  47     {
  48         label eI = pe[pointI][ppI];
  49
  50         const edge& e = es[eI];
  51
  52         if (e.start() == nextPointI || e.end() == nextPointI)
  53         {
  54             return eI;
  55         }
  56     }
  57
  58     // Make new edge.
  59     label edgeI = es.size();
  60     pe[pointI].append(edgeI);
  61     pe[nextPointI].append(edgeI);
  62     if (pointI < nextPointI)
  63     {
  64         es.append(edge(pointI, nextPointI));
  65     }
  66     else
  67     {
  68         es.append(edge(nextPointI, pointI));
  69     }
  70     return edgeI;
  71 }
  72
  73
  74 void Foam::primitiveMesh::calcEdges(const bool doFaceEdges) const
  75 {
  76     if (debug)
  77     {
  78         Pout<< "primitiveMesh::calcEdges(const bool) : "
  79             << "calculating edges, pointEdges and optionally faceEdges"
  80             << endl;
  81     }
  82
  83     // It is an error to attempt to recalculate edges
  84     // if the pointer is already set
  85     if ((edgesPtr_ || pePtr_) || (doFaceEdges && fePtr_))
  86     {
  87         FatalErrorIn("primitiveMesh::calcEdges(const bool) const")
  88             << "edges or pointEdges or faceEdges already calculated"
  89             << abort(FatalError);
  90     }
  91     else
  92     {
  93         // ALGORITHM:
  94         // Go through the faces list. Search pointEdges for existing edge.
  95         // If not found create edge and add to pointEdges.
  96         // In second loop sort edges in order of increasing neighbouring
  97         // point.
  98         // This algorithm replaces the one using pointFaces which used more
  99         // allocations but less memory and was on practical cases
 100         // quite a bit slower.
 101
 102         const faceList& fcs = faces();
 103
 104         // Size up lists
 105         // ~~~~~~~~~~~~~
 106
 107         // Estimate pointEdges storage
 108         List<DynamicList<label> > pe(nPoints());
 109         forAll(pe, pointI)
 110         {
 111             pe[pointI].setCapacity(primitiveMesh::edgesPerPoint_);
 112         }
 113
 114         // Estimate edges storage
 115         DynamicList<edge> es(pe.size()*primitiveMesh::edgesPerPoint_/2);
 116
 117         // Estimate faceEdges storage
 118         if (doFaceEdges)
 119         {
 120             fePtr_ = new labelListList(fcs.size());
 121             labelListList& faceEdges = *fePtr_;
 122             forAll(fcs, faceI)
 123             {
 124                 faceEdges[faceI].setSize(fcs[faceI].size());
 125             }
 126         }
 127
 128
 129         // Find consecutive face points in edge list
 130         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 131
 132         // Edges on boundary faces
 133         label nExtEdges = 0;
 134         // Edges using no boundary point
 135         nInternal0Edges_ = 0;
 136         // Edges using 1 boundary point
 137         label nInt1Edges = 0;
 138         // Edges using two boundary points but not on boundary face:
 139         // edges.size()-nExtEdges-nInternal0Edges_-nInt1Edges
 140
 141         // Ordering is different if points are ordered.
 142         if (nInternalPoints_ == -1)
 143         {
 144             // No ordering. No distinction between types.
 145             forAll(fcs, faceI)
 146             {
 147                 const face& f = fcs[faceI];
 148
 149                 forAll(f, fp)
 150                 {
 151                     label pointI = f[fp];
 152                     label nextPointI = f[f.fcIndex(fp)];
 153
 154                     label edgeI = getEdge(pe, es, pointI, nextPointI);
 155
 156                     if (doFaceEdges)
 157                     {
 158                         (*fePtr_)[faceI][fp] = edgeI;
 159                     }
 160                 }
 161             }
 162             // Assume all edges internal
 163             nExtEdges = 0;
 164             nInternal0Edges_ = es.size();
 165             nInt1Edges = 0;
 166         }
 167         else
 168         {
 169             // 1. Do external faces first. This creates external edges.
 170             for (label faceI = nInternalFaces_; faceI < fcs.size(); faceI++)
 171             {
 172                 const face& f = fcs[faceI];
 173
 174                 forAll(f, fp)
 175                 {
 176                     label pointI = f[fp];
 177                     label nextPointI = f[f.fcIndex(fp)];
 178
 179                     label oldNEdges = es.size();
 180                     label edgeI = getEdge(pe, es, pointI, nextPointI);
 181
 182                     if (es.size() > oldNEdges)
 183                     {
 184                         nExtEdges++;
 185                     }
 186                     if (doFaceEdges)
 187                     {
 188                         (*fePtr_)[faceI][fp] = edgeI;
 189                     }
 190                 }
 191             }
 192
 193             // 2. Do internal faces. This creates internal edges.
 194             for (label faceI = 0; faceI < nInternalFaces_; faceI++)
 195             {
 196                 const face& f = fcs[faceI];
 197
 198                 forAll(f, fp)
 199                 {
 200                     label pointI = f[fp];
 201                     label nextPointI = f[f.fcIndex(fp)];
 202
 203                     label oldNEdges = es.size();
 204                     label edgeI = getEdge(pe, es, pointI, nextPointI);
 205
 206                     if (es.size() > oldNEdges)
 207                     {
 208                         if (pointI < nInternalPoints_)
 209                         {
 210                             if (nextPointI < nInternalPoints_)
 211                             {
 212                                 nInternal0Edges_++;
 213                             }
 214                             else
 215                             {
 216                                 nInt1Edges++;
 217                             }
 218                         }
 219                         else
 220                         {
 221                             if (nextPointI < nInternalPoints_)
 222                             {
 223                                 nInt1Edges++;
 224                             }
 225                             else
 226                             {
 227                                 // Internal edge with two points on boundary
 228                             }
 229                         }
 230                     }
 231                     if (doFaceEdges)
 232                     {
 233                         (*fePtr_)[faceI][fp] = edgeI;
 234                     }
 235                 }
 236             }
 237         }
 238
 239
 240         // For unsorted meshes the edges will be in order of occurrence inside
 241         // the faces. For sorted meshes the first nExtEdges will be the external
 242         // edges.
 243
 244         if (nInternalPoints_ != -1)
 245         {
 246             nInternalEdges_ = es.size()-nExtEdges;
 247             nInternal1Edges_ = nInternal0Edges_+nInt1Edges;
 248
 249             //Pout<< "Edge overview:" << nl
 250             //    << "    total number of edges           : " << es.size()
 251             //    << nl
 252             //    << "    boundary edges                  : " << nExtEdges
 253             //    << nl
 254             //    << "    edges using no boundary point   : "
 255             //    << nInternal0Edges_
 256             //    << nl
 257             //    << "    edges using one boundary points : " << nInt1Edges
 258             //   << nl
 259             //    << "    edges using two boundary points : "
 260             //    << es.size()-nExtEdges-nInternal0Edges_-nInt1Edges << endl;
 261         }
 262
 263
 264         // Like faces sort edges in order of increasing neigbouring point.
 265         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 266         // Automatically if points are sorted into internal and external points
 267         // the edges will be sorted into
 268         // - internal point to internal point
 269         // - internal point to external point
 270         // - external point to external point
 271         // The problem is that the last one mixes external edges with internal
 272         // edges with two points on the boundary.
 273
 274
 275         // Map to sort into new upper-triangular order
 276         labelList oldToNew(es.size());
 277         if (debug)
 278         {
 279             oldToNew = -1;
 280         }
 281
 282         // start of edges with 0 boundary points
 283         label internal0EdgeI = 0;
 284
 285         // start of edges with 1 boundary points
 286         label internal1EdgeI = nInternal0Edges_;
 287
 288         // start of edges with 2 boundary points
 289         label internal2EdgeI = nInternal1Edges_;
 290
 291         // start of external edges
 292         label externalEdgeI = nInternalEdges_;
 293
 294
 295         // To sort neighbouring points in increasing order. Defined outside
 296         // for optimisation reasons: if all connectivity size same will need
 297         // no reallocations
 298         SortableList<label> nbrPoints(primitiveMesh::edgesPerPoint_);
 299
 300         forAll(pe, pointI)
 301         {
 302             const DynamicList<label>& pEdges = pe[pointI];
 303
 304             nbrPoints.setSize(pEdges.size());
 305
 306             forAll(pEdges, i)
 307             {
 308                 const edge& e = es[pEdges[i]];
 309
 310                 label nbrPointI = e.otherVertex(pointI);
 311
 312                 if (nbrPointI < pointI)
 313                 {
 314                     nbrPoints[i] = -1;
 315                 }
 316                 else
 317                 {
 318                     nbrPoints[i] = nbrPointI;
 319                 }
 320             }
 321             nbrPoints.sort();
 322
 323
 324             if (nInternalPoints_ == -1)
 325             {
 326                 // Sort all upper-triangular
 327                 forAll(nbrPoints, i)
 328                 {
 329                     if (nbrPoints[i] != -1)
 330                     {
 331                         label edgeI = pEdges[nbrPoints.indices()[i]];
 332
 333                         oldToNew[edgeI] = internal0EdgeI++;
 334                     }
 335                 }
 336             }
 337             else
 338             {
 339                 if (pointI < nInternalPoints_)
 340                 {
 341                     forAll(nbrPoints, i)
 342                     {
 343                         label nbrPointI = nbrPoints[i];
 344
 345                         label edgeI = pEdges[nbrPoints.indices()[i]];
 346
 347                         if (nbrPointI != -1)
 348                         {
 349                             if (edgeI < nExtEdges)
 350                             {
 351                                 // External edge
 352                                 oldToNew[edgeI] = externalEdgeI++;
 353                             }
 354                             else if (nbrPointI < nInternalPoints_)
 355                             {
 356                                 // Both points inside
 357                                 oldToNew[edgeI] = internal0EdgeI++;
 358                             }
 359                             else
 360                             {
 361                                 // One points inside, one outside
 362                                 oldToNew[edgeI] = internal1EdgeI++;
 363                             }
 364                         }
 365                     }
 366                 }
 367                 else
 368                 {
 369                     forAll(nbrPoints, i)
 370                     {
 371                         label nbrPointI = nbrPoints[i];
 372
 373                         label edgeI = pEdges[nbrPoints.indices()[i]];
 374
 375                         if (nbrPointI != -1)
 376                         {
 377                             if (edgeI < nExtEdges)
 378                             {
 379                                 // External edge
 380                                 oldToNew[edgeI] = externalEdgeI++;
 381                             }
 382                             else if (nbrPointI < nInternalPoints_)
 383                             {
 384                                 // Not possible!
 385                                 FatalErrorIn("primitiveMesh::calcEdges(..)")
 386                                     << abort(FatalError);
 387                             }
 388                             else
 389                             {
 390                                 // Both points outside
 391                                 oldToNew[edgeI] = internal2EdgeI++;
 392                             }
 393                         }
 394                     }
 395                 }
 396             }
 397         }
 398
 399
 400         if (debug)
 401         {
 402             label edgeI = findIndex(oldToNew, -1);
 403
 404             if (edgeI != -1)
 405             {
 406                 const edge& e = es[edgeI];
 407
 408                 FatalErrorIn("primitiveMesh::calcEdges(..)")
 409                     << "Did not sort edge " << edgeI << " points:" << e
 410                     << " coords:" << points()[e[0]] << points()[e[1]]
 411                     << endl
 412                     << "Current buckets:" << endl
 413                     << "    internal0EdgeI:" << internal0EdgeI << endl
 414                     << "    internal1EdgeI:" << internal1EdgeI << endl
 415                     << "    internal2EdgeI:" << internal2EdgeI << endl
 416                     << "    externalEdgeI:" << externalEdgeI << endl
 417                     << abort(FatalError);
 418             }
 419         }
 420
 421
 422
 423         // Renumber and transfer.
 424
 425         // Edges
 426         edgesPtr_ = new edgeList(es.size());
 427         edgeList& edges = *edgesPtr_;
 428         forAll(es, edgeI)
 429         {
 430             edges[oldToNew[edgeI]] = es[edgeI];
 431         }
 432
 433         // pointEdges
 434         pePtr_ = new labelListList(nPoints());
 435         labelListList& pointEdges = *pePtr_;
 436         forAll(pe, pointI)
 437         {
 438             DynamicList<label>& pEdges = pe[pointI];
 439             pEdges.shrink();
 440             inplaceRenumber(oldToNew, pEdges);
 441             pointEdges[pointI].transfer(pEdges);
 442             Foam::sort(pointEdges[pointI]);
 443         }
 444
 445         // faceEdges
 446         if (doFaceEdges)
 447         {
 448             labelListList& faceEdges = *fePtr_;
 449             forAll(faceEdges, faceI)
 450             {
 451                 inplaceRenumber(oldToNew, faceEdges[faceI]);
 452             }
 453         }
 454     }
 455 }
 456
 457
 458 Foam::label Foam::primitiveMesh::findFirstCommonElementFromSortedLists
 459 (
 460     const labelList& list1,
 461     const labelList& list2
 462 )
 463 {
 464     label result = -1;
 465
 466     labelList::const_iterator iter1 = list1.begin();
 467     labelList::const_iterator iter2 = list2.begin();
 468
 469     while (iter1 != list1.end() && iter2 != list2.end())
 470     {
 471         if( *iter1 < *iter2)
 472         {
 473             ++iter1;
 474         }
 475         else if (*iter1 > *iter2)
 476         {
 477             ++iter2;
 478         }
 479         else
 480         {
 481             result = *iter1;
 482             break;
 483         }
 484     }
 485     if (result == -1)
 486     {
 487         FatalErrorIn
 488         (
 489             "primitiveMesh::findFirstCommonElementFromSortedLists"
 490             "(const labelList&, const labelList&)"
 491         )   << "No common elements in lists " << list1 << " and " << list2
 492             << abort(FatalError);
 493     }
 494     return result;
 495 }
 496
 497
 498 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 499
 500 const Foam::edgeList& Foam::primitiveMesh::edges() const
 501 {
 502     if (!edgesPtr_)
 503     {
 504         //calcEdges(true);
 505         calcEdges(false);
 506     }
 507
 508     return *edgesPtr_;
 509 }
 510
 511 const Foam::labelListList& Foam::primitiveMesh::pointEdges() const
 512 {
 513     if (!pePtr_)
 514     {
 515         //calcEdges(true);
 516         calcEdges(false);
 517     }
 518
 519     return *pePtr_;
 520 }
 521
 522
 523 const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
 524 {
 525     if (!fePtr_)
 526     {
 527         if (debug)
 528         {
 529             Pout<< "primitiveMesh::faceEdges() : "
 530                 << "calculating faceEdges" << endl;
 531         }
 532
 533         //calcEdges(true);
 534         const faceList& fcs = faces();
 535         const labelListList& pe = pointEdges();
 536         const edgeList& es = edges();
 537
 538         fePtr_ = new labelListList(fcs.size());
 539         labelListList& faceEdges = *fePtr_;
 540
 541         forAll(fcs, faceI)
 542         {
 543             const face& f = fcs[faceI];
 544
 545             labelList& fEdges = faceEdges[faceI];
 546             fEdges.setSize(f.size());
 547
 548             forAll(f, fp)
 549             {
 550                 label pointI = f[fp];
 551                 label nextPointI = f[f.fcIndex(fp)];
 552
 553                 // Find edge between pointI, nextPontI
 554                 const labelList& pEdges = pe[pointI];
 555
 556                 forAll(pEdges, i)
 557                 {
 558                     label edgeI = pEdges[i];
 559
 560                     if (es[edgeI].otherVertex(pointI) == nextPointI)
 561                     {
 562                         fEdges[fp] = edgeI;
 563                         break;
 564                     }
 565                 }
 566             }
 567         }
 568     }
 569
 570     return *fePtr_;
 571 }
 572
 573
 574 void Foam::primitiveMesh::clearOutEdges()
 575 {
 576     deleteDemandDrivenData(edgesPtr_);
 577     deleteDemandDrivenData(pePtr_);
 578     deleteDemandDrivenData(fePtr_);
 579     labels_.clear();
 580     labelSet_.clear();
 581 }
 582
 583
 584 const Foam::labelList& Foam::primitiveMesh::faceEdges
 585 (
 586     const label faceI,
 587     DynamicList<label>& storage
 588 ) const
 589 {
 590     if (hasFaceEdges())
 591     {
 592         return faceEdges()[faceI];
 593     }
 594     else
 595     {
 596         const labelListList& pointEs = pointEdges();
 597         const face& f = faces()[faceI];
 598
 599         storage.clear();
 600         if (f.size() > storage.capacity())
 601         {
 602             storage.setCapacity(f.size());
 603         }
 604
 605         forAll(f, fp)
 606         {
 607             storage.append
 608             (
 609                 findFirstCommonElementFromSortedLists
 610                 (
 611                     pointEs[f[fp]],
 612                     pointEs[f.nextLabel(fp)]
 613                 )
 614             );
 615         }
 616
 617         return storage;
 618     }
 619 }
 620
 621
 622 const Foam::labelList& Foam::primitiveMesh::faceEdges(const label faceI) const
 623 {
 624     return faceEdges(faceI, labels_);
 625 }
 626
 627
 628 const Foam::labelList& Foam::primitiveMesh::cellEdges
 629 (
 630     const label cellI,
 631     DynamicList<label>& storage
 632 ) const
 633 {
 634     if (hasCellEdges())
 635     {
 636         return cellEdges()[cellI];
 637     }
 638     else
 639     {
 640         const labelList& cFaces = cells()[cellI];
 641
 642         labelSet_.clear();
 643
 644         forAll(cFaces, i)
 645         {
 646             const labelList& fe = faceEdges(cFaces[i]);
 647
 648             forAll(fe, feI)
 649             {
 650                 labelSet_.insert(fe[feI]);
 651             }
 652         }
 653
 654         storage.clear();
 655
 656         if (labelSet_.size() > storage.capacity())
 657         {
 658             storage.setCapacity(labelSet_.size());
 659         }
 660
 661         forAllConstIter(labelHashSet, labelSet_, iter)
 662         {
 663             storage.append(iter.key());
 664         }
 665
 666         return storage;
 667     }
 668 }
 669
 670
 671 const Foam::labelList& Foam::primitiveMesh::cellEdges(const label cellI) const
 672 {
 673     return cellEdges(cellI, labels_);
 674 }
 675
 676
 677 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 678
 679 // ************************************************************************* //