src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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
  75 {
  76     if (debug)
  77     {
  78         Pout<< "primitiveMesh::calcEdges() : "
  79             << "calculating edges, pointEdges and 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_ || fePtr_)
  86     {
  87         FatalErrorIn("primitiveMesh::calcEdges(const bool) const")
  88             << "edges or faceEdges already calculated"
  89             << abort(FatalError);
  90     }
  91     else
  92     {
  93         // Note: replaced with the original and correct algorithm
  94         // which preserved correct ordering for edges list
  95         // Version 1.5.x, 1.5-dev were wrong, with failures in parallel
  96         // point-based data exchange.  Fixed in 1.6-ext and
  97         // consecutive versions
  98         // HJ, 27/Aug/2010
  99
 100         // ALGORITHM:
 101         // Go through the pointFace list.  Go through the list of faces for that
 102         // point and ask for edges.  If the edge has got the point in question
 103         // AND the second point in the edge is larger than the first, add the
 104         // edge to the list.  At the same time, add the edge label to the list
 105         // of edges for the current face (faceEdges) and log the other face as
 106         // the neighbour of this face.
 107
 108         const faceList& f = faces();
 109
 110         const labelListList& pf = pointFaces();
 111
 112         fePtr_ = new labelListList(nFaces());
 113         labelListList& fe = *fePtr_;
 114
 115         // count the maximum number of edges
 116         label maxEdges = 0;
 117
 118         // create a list of edges for each face and store for efficiency
 119         edgeListList edgesOfFace(nFaces());
 120
 121         forAll (f, faceI)
 122         {
 123             edgesOfFace[faceI] = f[faceI].edges();
 124
 125             maxEdges += f[faceI].nEdges();
 126
 127             labelList& curFE = fe[faceI];
 128
 129             curFE.setSize(f[faceI].nEdges());
 130
 131             forAll (curFE, curFEI)
 132             {
 133                 curFE[curFEI] = -1;
 134             }
 135         }
 136
 137         // EDGE CALCULATION
 138
 139         edgesPtr_ = new edgeList(maxEdges);
 140         edgeList& e = *edgesPtr_;
 141         label nEdges = 0;
 142
 143         forAll (pf, pointI)
 144         {
 145             const labelList& curFaces = pf[pointI];
 146
 147             // create a list of labels to keep the neighbours that
 148             // have already been added
 149             DynamicList<label, edgesPerPoint_> addedNeighbours;
 150             DynamicList<DynamicList<label, edgesPerPoint_> >
 151                 faceGivingNeighbour;
 152             DynamicList<DynamicList<label, edgesPerPoint_> >
 153                 edgeOfFaceGivingNeighbour;
 154
 155             forAll (curFaces, faceI)
 156             {
 157                 // get the edges
 158                 const edgeList& fEdges = edgesOfFace[curFaces[faceI]];
 159
 160                 // for every edge
 161                 forAll(fEdges, edgeI)
 162                 {
 163                     const edge& ends = fEdges[edgeI];
 164
 165                     // does the edge has got the point in question
 166                     bool found = false;
 167                     label secondPoint = -1;
 168
 169                     if (ends.start() == pointI)
 170                     {
 171                         found = true;
 172                         secondPoint = ends.end();
 173                     }
 174
 175                     if (ends.end() == pointI)
 176                     {
 177                         found = true;
 178                         secondPoint = ends.start();
 179                     }
 180
 181                     // if the edge has got the point and second label is larger
 182                     // than first, it is a candidate for adding
 183                     if (found && (secondPoint > pointI))
 184                     {
 185                         // check if the edge has already been added
 186                         bool added = false;
 187
 188                         forAll (addedNeighbours, eopI)
 189                         {
 190                             if (secondPoint == addedNeighbours[eopI])
 191                             {
 192                                 // Edge is already added. New face sharing it
 193                                 added = true;
 194
 195                                 // Remember the face and edge giving neighbour
 196                                 faceGivingNeighbour[eopI].append
 197                                     (curFaces[faceI]);
 198
 199                                 edgeOfFaceGivingNeighbour[eopI].append(edgeI);
 200
 201                                 break;
 202                             }
 203                         }
 204
 205                         // If not added, add the edge to the list
 206                         if (!added)
 207                         {
 208                             addedNeighbours.append(secondPoint);
 209
 210                             // Remember the face and subShape giving neighbour
 211                             faceGivingNeighbour(addedNeighbours.size() - 1)
 212                                 .append(curFaces[faceI]);
 213                             edgeOfFaceGivingNeighbour
 214                                 (addedNeighbours.size() - 1).append(edgeI);
 215                         }
 216                     }
 217                 }
 218             }
 219
 220             // All edges for the current point found. Before adding them to the
 221             // list, it is necessary to sort them in the increasing order of the
 222             // neighbouring point.
 223
 224             // Make real list out of SLList to simplify the manipulation.
 225             // Also, make another list to "remember" how the original list was
 226             // reshuffled.
 227             labelList shuffleList(addedNeighbours.size());
 228
 229             forAll (shuffleList, i)
 230             {
 231                 shuffleList[i] = i;
 232             }
 233
 234             // Use a simple sort to sort the addedNeighbours list.
 235             //  Other two lists mimic the same behaviour
 236             label i, j, a, b;
 237
 238             for (j = 1; j <= addedNeighbours.size() - 1; j++)
 239             {
 240                 a = addedNeighbours[j];
 241                 b = shuffleList[j];
 242
 243                 i = j - 1;
 244
 245                 while (i >= 0 && addedNeighbours[i] > a)
 246                 {
 247                     addedNeighbours[i + 1] = addedNeighbours[i];
 248                     shuffleList[i + 1] = shuffleList[i];
 249                     i--;
 250                 }
 251
 252                 addedNeighbours[i + 1] = a;
 253                 shuffleList[i + 1] = b;
 254             }
 255
 256             labelList reshuffleList(shuffleList.size());
 257
 258             forAll(shuffleList, i)
 259             {
 260                 reshuffleList[shuffleList[i]] = i;
 261             }
 262
 263             // Reshuffle other lists
 264
 265             labelListList fgn(faceGivingNeighbour.size());
 266
 267             forAll (faceGivingNeighbour, i)
 268             {
 269                 fgn[reshuffleList[i]].transfer(faceGivingNeighbour[i].shrink());
 270             }
 271
 272             labelListList eofgn(edgeOfFaceGivingNeighbour.size());
 273
 274             forAll (edgeOfFaceGivingNeighbour, i)
 275             {
 276                 eofgn[reshuffleList[i]].transfer
 277                 (
 278                     edgeOfFaceGivingNeighbour[i].shrink()
 279                 );
 280             }
 281
 282             // adding the edges
 283             forAll(addedNeighbours, edgeI)
 284             {
 285                 const labelList& curFgn = fgn[edgeI];
 286                 const labelList& curEofgn = eofgn[edgeI];
 287
 288                 forAll (curFgn, fgnI)
 289                 {
 290                     fe[curFgn[fgnI]][curEofgn[fgnI]] = nEdges;
 291                 }
 292
 293                 e[nEdges] = edge(pointI, addedNeighbours[edgeI]);
 294                 nEdges++;
 295             }
 296         }
 297
 298         // reset the size
 299         e.setSize(nEdges);
 300     }
 301 }
 302
 303 // Removed ordered edges: reverting to correct parallelisation
 304 // of edge data.  HJ, 17/Au/2010
 305 // void primitiveMesh::calcOrderedEdges() const
 306
 307
 308 Foam::label Foam::primitiveMesh::findFirstCommonElementFromSortedLists
 309 (
 310     const labelList& list1,
 311     const labelList& list2
 312 )
 313 {
 314     label result = -1;
 315
 316     labelList::const_iterator iter1 = list1.begin();
 317     labelList::const_iterator iter2 = list2.begin();
 318
 319     while (iter1 != list1.end() && iter2 != list2.end())
 320     {
 321         if( *iter1 < *iter2)
 322         {
 323             ++iter1;
 324         }
 325         else if (*iter1 > *iter2)
 326         {
 327             ++iter2;
 328         }
 329         else
 330         {
 331             result = *iter1;
 332             break;
 333         }
 334     }
 335     if (result == -1)
 336     {
 337         FatalErrorIn
 338         (
 339             "primitiveMesh::findFirstCommonElementFromSortedLists"
 340             "(const labelList&, const labelList&)"
 341         )   << "No common elements in lists " << list1 << " and " << list2
 342             << abort(FatalError);
 343     }
 344     return result;
 345 }
 346
 347
 348 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 349
 350 const Foam::edgeList& Foam::primitiveMesh::edges() const
 351 {
 352     if (!edgesPtr_)
 353     {
 354         // 1.6.x merge: reverting to correct code
 355         // HJ, 17/Aug/2010
 356         calcEdges();
 357     }
 358
 359     return *edgesPtr_;
 360 }
 361
 362
 363 const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
 364 {
 365     if (!fePtr_)
 366     {
 367         // 1.6.x merge: reverting to correct code
 368         // HJ, 17/Aug/2010
 369         calcEdges();
 370     }
 371
 372     return *fePtr_;
 373 }
 374
 375
 376 void Foam::primitiveMesh::clearOutEdges()
 377 {
 378     deleteDemandDrivenData(edgesPtr_);
 379     deleteDemandDrivenData(fePtr_);
 380     labels_.clear();
 381     labelSet_.clear();
 382 }
 383
 384
 385 const Foam::labelList& Foam::primitiveMesh::faceEdges
 386 (
 387     const label faceI,
 388     DynamicList<label>& storage
 389 ) const
 390 {
 391     if (hasFaceEdges())
 392     {
 393         return faceEdges()[faceI];
 394     }
 395     else
 396     {
 397         const labelListList& pointEs = pointEdges();
 398         const face& f = faces()[faceI];
 399
 400         storage.clear();
 401         if (f.size() > storage.capacity())
 402         {
 403             storage.setCapacity(f.size());
 404         }
 405
 406         forAll(f, fp)
 407         {
 408             storage.append
 409             (
 410                 findFirstCommonElementFromSortedLists
 411                 (
 412                     pointEs[f[fp]],
 413                     pointEs[f.nextLabel(fp)]
 414                 )
 415             );
 416         }
 417
 418         return storage;
 419     }
 420 }
 421
 422
 423 const Foam::labelList& Foam::primitiveMesh::faceEdges(const label faceI) const
 424 {
 425     return faceEdges(faceI, labels_);
 426 }
 427
 428
 429 const Foam::labelList& Foam::primitiveMesh::cellEdges
 430 (
 431     const label cellI,
 432     DynamicList<label>& storage
 433 ) const
 434 {
 435     if (hasCellEdges())
 436     {
 437         return cellEdges()[cellI];
 438     }
 439     else
 440     {
 441         const labelList& cFaces = cells()[cellI];
 442
 443         labelSet_.clear();
 444
 445         forAll(cFaces, i)
 446         {
 447             const labelList& fe = faceEdges(cFaces[i]);
 448
 449             forAll(fe, feI)
 450             {
 451                 labelSet_.insert(fe[feI]);
 452             }
 453         }
 454
 455         storage.clear();
 456
 457         if (labelSet_.size() > storage.capacity())
 458         {
 459             storage.setCapacity(labelSet_.size());
 460         }
 461
 462         forAllConstIter(labelHashSet, labelSet_, iter)
 463         {
 464             storage.append(iter.key());
 465         }
 466
 467         return storage;
 468     }
 469 }
 470
 471
 472 const Foam::labelList& Foam::primitiveMesh::cellEdges(const label cellI) const
 473 {
 474     return cellEdges(cellI, labels_);
 475 }
 476
 477
 478 // ************************************************************************* //