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