src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAddressing/meshOctreeAddressingCreation.C

   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 "meshOctreeAddressing.H"
  29 #include "helperFunctions.H"
  30 #include "VRWGraphSMPModifier.H"
  31 #include "demandDrivenData.H"
  32 #include "meshOctree.H"
  33 #include "labelLongList.H"
  34 #include "triSurf.H"
  35
  36 # ifdef USE_OMP
  37 #include <omp.h>
  38 # endif
  39
  40 //#define DEBUGVrt
  41
  42 # ifdef DEBUGVrt
  43 #include "OFstream.H"
  44 # endif
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 namespace Foam
  49 {
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53 # ifdef DEBUGVrt
  54 void writeOctreeToVTK
  55 (
  56     const fileName& fName,
  57     const meshOctree& octree,
  58     const List<direction>& boxTypes,
  59     const direction bType
  60 )
  61 {
  62     OFstream file(fName);
  63
  64     //- write the header
  65     file << "# vtk DataFile Version 3.0\n";
  66     file << "vtk output\n";
  67     file << "ASCII\n";
  68     file << "DATASET UNSTRUCTURED_GRID\n";
  69
  70     label nBoxes(0);
  71     forAll(boxTypes, leafI)
  72         if( boxTypes[leafI] & bType )
  73             ++nBoxes;
  74
  75     //- write points
  76     file << "POINTS " << 8*nBoxes << " float\n";
  77     forAll(boxTypes, leafI)
  78     {
  79         if( boxTypes[leafI] & bType )
  80         {
  81             FixedList<point, 8> vertices;
  82             octree.returnLeaf(leafI).vertices(octree.rootBox(), vertices);
  83
  84             forAll(vertices, vI)
  85             {
  86                 const point& p = vertices[vI];
  87
  88                 file << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
  89             }
  90         }
  91     }
  92
  93     //- write boxes
  94     file << "\nCELLS " << nBoxes
  95          << " " << 9 * nBoxes << nl;
  96
  97     nBoxes = 0;
  98     forAll(boxTypes, leafI)
  99     {
 100         if( boxTypes[leafI] & bType )
 101         {
 102             const label start = 8 * nBoxes;
 103             file << 8 << " " << start << " " << start+1
 104                       << " " << start+3 << " " << start+2
 105                       << " " << start+4 << " " << start+5
 106                       << " " << start+7 << " " << start+6 << nl;
 107
 108             ++nBoxes;
 109         }
 110     }
 111
 112     file << nl;
 113
 114     //- write cell types
 115     file << "CELL_TYPES " << nBoxes << nl;
 116     for(label i=0;i<nBoxes;++i)
 117         file << 12 << nl;
 118
 119     file << nl;
 120 }
 121 # endif
 122
 123 void meshOctreeAddressing::createOctreePoints() const
 124 {
 125     const VRWGraph& nodeLabels = this->nodeLabels();
 126     const boundBox& rootBox = octree_.rootBox();
 127
 128     octreePointsPtr_ = new pointField(nNodes_);
 129     pointField& octreePoints = *octreePointsPtr_;
 130
 131     const label nLeaves = nodeLabels.size();
 132     # ifdef USE_OMP
 133     # pragma omp parallel for schedule(guided, 100)
 134     # endif
 135     for(label cubeI=0;cubeI<nLeaves;++cubeI)
 136     {
 137         if( nodeLabels.sizeOfRow(cubeI) == 0 )
 138             continue;
 139
 140         FixedList<point, 8> vertices;
 141         const meshOctreeCubeBasic& oc = octree_.returnLeaf(cubeI);
 142         oc.vertices(rootBox, vertices);
 143
 144         forAllRow(nodeLabels, cubeI, nI)
 145         {
 146             const label nodeI = nodeLabels(cubeI, nI);
 147
 148             octreePoints[nodeI] = vertices[nI];
 149         }
 150     }
 151 }
 152
 153 void meshOctreeAddressing::createNodeLabels() const
 154 {
 155     const List<direction>& boxType = this->boxType();
 156
 157     nodeLabelsPtr_ = new VRWGraph(octree_.numberOfLeaves());
 158     VRWGraph& nodeLabels = *nodeLabelsPtr_;
 159
 160     //- allocate storage for node labels
 161     forAll(nodeLabels, leafI)
 162     {
 163         if( boxType[leafI] )
 164         {
 165             nodeLabels.setRowSize(leafI, 8);
 166
 167             forAllRow(nodeLabels, leafI, i)
 168                 nodeLabels(leafI, i) = -1;
 169         }
 170     }
 171
 172     //- start creating node labels
 173     nNodes_ = 0;
 174     DynList<label> numLocalNodes;
 175     # ifdef USE_OMP
 176     # pragma omp parallel
 177     # endif
 178     {
 179         # ifdef USE_OMP
 180         const label nThreads = omp_get_num_threads();
 181         const label threadI = omp_get_thread_num();
 182         # else
 183         const label nThreads = 1;
 184         const label threadI = 0;
 185         # endif
 186
 187         # ifdef USE_OMP
 188         # pragma omp master
 189         # endif
 190         numLocalNodes.setSize(nThreads);
 191
 192         # ifdef USE_OMP
 193         # pragma omp barrier
 194         # endif
 195
 196         //- count the number of nodes local to each process
 197         label& nLocalNodes = numLocalNodes[threadI];
 198         nLocalNodes = 0;
 199
 200         # ifdef USE_OMP
 201         # pragma omp for schedule(static, 100)
 202         # endif
 203         forAll(nodeLabels, leafI)
 204         {
 205             forAllRow(nodeLabels, leafI, nI)
 206             {
 207                 if( nodeLabels(leafI, nI) != -1 )
 208                     continue;
 209
 210                 FixedList<label, 8> pLeaves;
 211                 octree_.findLeavesForCubeVertex(leafI, nI, pLeaves);
 212
 213                 FixedList<bool, 8> validLeaf(true);
 214                 label minLeaf(leafI);
 215                 forAll(pLeaves, plI)
 216                 {
 217                     if( pLeaves[plI] > -1 )
 218                     {
 219                         for(label i=plI+1;i<8;++i)
 220                             if( pLeaves[plI] == pLeaves[i] )
 221                             {
 222                                 validLeaf[plI] = false;
 223                                 validLeaf[i] = false;
 224                             }
 225
 226                         if( !boxType[pLeaves[plI]] )
 227                         {
 228                             validLeaf[plI] = false;
 229                             pLeaves[plI] = -1;
 230                         }
 231
 232                         if( validLeaf[plI] )
 233                             minLeaf = Foam::min(minLeaf, pLeaves[plI]);
 234                     }
 235                     else
 236                     {
 237                         validLeaf[plI] = false;
 238                     }
 239                 }
 240
 241                 if( (minLeaf == leafI) && validLeaf[7-nI] )
 242                 {
 243                     forAll(pLeaves, plI)
 244                         if( validLeaf[plI] )
 245                         {
 246                             //- set node labels to -2 not to repeat searches
 247                             nodeLabels(pLeaves[plI], (7-plI)) = -2;
 248                         }
 249
 250                     ++nLocalNodes;
 251                 }
 252             }
 253         }
 254
 255         //- set start node for each process
 256         # ifdef USE_OMP
 257         # pragma omp barrier
 258         # endif
 259
 260         label startNode(0);
 261         for(label i=0;i<threadI;++i)
 262             startNode += numLocalNodes[i];
 263
 264         //- start creating node labels
 265         # ifdef USE_OMP
 266         # pragma omp for schedule(static, 100)
 267         # endif
 268         forAll(nodeLabels, leafI)
 269         {
 270             forAllRow(nodeLabels, leafI, nI)
 271             {
 272                 if( nodeLabels(leafI, nI) >= 0 )
 273                     continue;
 274
 275                 FixedList<label, 8> pLeaves;
 276                 octree_.findLeavesForCubeVertex(leafI, nI, pLeaves);
 277
 278                 FixedList<bool, 8> validLeaf(true);
 279                 label minLeaf(leafI);
 280                 forAll(pLeaves, plI)
 281                 {
 282                     if( pLeaves[plI] > -1 )
 283                     {
 284                         for(label i=plI+1;i<8;++i)
 285                             if( pLeaves[plI] == pLeaves[i] )
 286                             {
 287                                 validLeaf[plI] = false;
 288                                 validLeaf[i] = false;
 289                             }
 290
 291                         if( !boxType[pLeaves[plI]] )
 292                         {
 293                             validLeaf[plI] = false;
 294                             pLeaves[plI] = -1;
 295                         }
 296
 297                         if( validLeaf[plI] )
 298                             minLeaf = Foam::min(minLeaf, pLeaves[plI]);
 299                     }
 300                     else
 301                     {
 302                         validLeaf[plI] = false;
 303                     }
 304                 }
 305
 306                 if( (minLeaf == leafI) && validLeaf[7-nI] )
 307                 {
 308                     forAll(pLeaves, plI)
 309                         if( validLeaf[plI] )
 310                         {
 311                             //- store the vertex at the corresponding
 312                             //- location in the cube
 313                             nodeLabels(pLeaves[plI], (7-plI)) = startNode;
 314                         }
 315
 316                     //- store vertex label
 317                     ++startNode;
 318                 }
 319             }
 320         }
 321
 322         //- set the number of nodes
 323         # ifdef USE_OMP
 324         # pragma omp critical
 325         # endif
 326         {
 327             nNodes_ = Foam::max(nNodes_, startNode);
 328         }
 329     }
 330
 331     # ifdef DEBUGVrt
 332     List<direction> badLeaves(nodeLabels.size(), direction(0));
 333     forAll(nodeLabels, leafI)
 334         forAllRow(nodeLabels, leafI, i)
 335             if( nodeLabels(leafI, i) < 0 )
 336                 badLeaves[leafI] |= 1;
 337     writeOctreeToVTK("badLeaves.vtk", octree_, badLeaves, 1);
 338
 339     writeOctreeToVTK("meshCells.vtk", octree_, boxType, MESHCELL);
 340     writeOctreeToVTK("boundaryCells.vtk", octree_, boxType, BOUNDARY);
 341
 342     Info << "Checking for existence of negative node labels" << endl;
 343     forAll(nodeLabels, leafI)
 344     {
 345         forAllRow(nodeLabels, leafI, nI)
 346         {
 347             if( nodeLabels(leafI, nI) < 0 )
 348             {
 349                 FixedList<label, 8> pLeaves;
 350                 octree_.findLeavesForCubeVertex(leafI, nI, pLeaves);
 351
 352                 FixedList<bool, 8> validLeaf(true);
 353                 label minLeaf(leafI);
 354                 forAll(pLeaves, plI)
 355                 {
 356                     if( pLeaves[plI] > -1 )
 357                     {
 358                         for(label i=plI+1;i<8;++i)
 359                             if( pLeaves[plI] == pLeaves[i] )
 360                             {
 361                                 validLeaf[plI] = false;
 362                                 validLeaf[i] = false;
 363                             }
 364
 365                         if( !boxType[pLeaves[plI]] )
 366                         {
 367                             validLeaf[plI] = false;
 368                             pLeaves[plI] = -1;
 369                         }
 370
 371                         if( validLeaf[plI] )
 372                             minLeaf = Foam::min(minLeaf, pLeaves[plI]);
 373                     }
 374                     else
 375                     {
 376                         validLeaf[plI] = false;
 377                     }
 378                 }
 379
 380                 Info << "Min leaf " << minLeaf << endl;
 381                 Info << "Valid leaf " << validLeaf << endl;
 382                 Info << "pLeaves " << pLeaves << endl;
 383                 Info << "Node position " << nI << endl;
 384
 385                 Info << "1.Leaf " << leafI << " node labels "
 386                       << nodeLabels[leafI] << endl;
 387
 388                 forAll(validLeaf, i)
 389                     if( validLeaf[i] )
 390                         Info << "Leaf at position " << i << " has node labels "
 391                              << nodeLabels[pLeaves[i]]
 392                              << " at level "
 393                              << octree_.returnLeaf(pLeaves[i]).level() << endl;
 394             }
 395         }
 396     }
 397     # endif
 398 }
 399
 400 void meshOctreeAddressing::createNodeLeaves() const
 401 {
 402     const List<direction>& boxType = this->boxType();
 403     const VRWGraph& nodeLabels = this->nodeLabels();
 404
 405     //- allocate nodeLeavesPtr_
 406     nodeLeavesPtr_ = new FRWGraph<label, 8>(nNodes_);
 407     FRWGraph<label, 8>& nodeLeaves = *nodeLeavesPtr_;
 408
 409     boolList storedNode(nNodes_, false);
 410     # ifdef USE_OMP
 411     # pragma omp parallel for schedule(dynamic, 100)
 412     # endif
 413     forAll(nodeLabels, leafI)
 414     {
 415         forAllRow(nodeLabels, leafI, nI)
 416         {
 417             const label nodeI = nodeLabels(leafI, nI);
 418
 419             if( storedNode[nodeI] )
 420                 continue;
 421
 422             storedNode[nodeI] = true;
 423
 424             FixedList<label, 8> pLeaves;
 425             octree_.findLeavesForCubeVertex(leafI, nI, pLeaves);
 426
 427             forAll(pLeaves, plI)
 428             {
 429                 if( pLeaves[plI] < 0 )
 430                     continue;
 431
 432                 if( !boxType[pLeaves[plI]] )
 433                     pLeaves[plI] = -1;
 434             }
 435
 436             nodeLeaves.setRow(nodeI, pLeaves);
 437         }
 438     }
 439 }
 440
 441 void meshOctreeAddressing::findUsedBoxes() const
 442 {
 443     boxTypePtr_ = new List<direction>(octree_.numberOfLeaves(), NONE);
 444     List<direction>& boxType = *boxTypePtr_;
 445
 446     # ifdef USE_OMP
 447     # pragma omp parallel for schedule(dynamic, 40)
 448     # endif
 449     forAll(boxType, leafI)
 450     {
 451         const meshOctreeCubeBasic& leaf = octree_.returnLeaf(leafI);
 452
 453         if(
 454             !octree_.hasContainedTriangles(leafI) &&
 455             !octree_.hasContainedEdges(leafI) &&
 456             (leaf.cubeType() & meshOctreeCubeBasic::INSIDE)
 457         )
 458             boxType[leafI] |= MESHCELL;
 459     }
 460
 461     if( meshDict_.found("nonManifoldMeshing") )
 462     {
 463         const bool nonManifoldMesh
 464         (
 465             readBool(meshDict_.lookup("nonManifoldMeshing"))
 466         );
 467
 468         if( nonManifoldMesh )
 469         {
 470             # ifdef USE_OMP
 471             # pragma omp parallel for schedule(dynamic, 40)
 472             # endif
 473             forAll(boxType, leafI)
 474             {
 475                 const meshOctreeCubeBasic& leaf = octree_.returnLeaf(leafI);
 476                     if( leaf.cubeType() & meshOctreeCubeBasic::UNKNOWN )
 477                         boxType[leafI] |= MESHCELL;
 478             }
 479         }
 480     }
 481
 482     if( useDATABoxes_ )
 483     {
 484         Info << "Using DATA boxes" << endl;
 485
 486         forAll(boxType, leafI)
 487         {
 488             if(
 489                 octree_.hasContainedTriangles(leafI) ||
 490                 octree_.hasContainedEdges(leafI)
 491             )
 492                 boxType[leafI] |= MESHCELL;
 493         }
 494
 495         //- do not use boxes intersecting given patches
 496         if( meshDict_.found("removeCellsIntersectingPatches") )
 497         {
 498             wordHashSet patchesToRemove;
 499
 500             if( meshDict_.isDict("removeCellsIntersectingPatches") )
 501             {
 502                 const dictionary& dict =
 503                     meshDict_.subDict("removeCellsIntersectingPatches");
 504                 const wordList patchNames = dict.toc();
 505                 forAll(patchNames, patchI)
 506                     patchesToRemove.insert(patchNames[patchI]);
 507             }
 508             else
 509             {
 510                 wordHashSet patchesToRemoveCopy
 511                 (
 512                     meshDict_.lookup("removeCellsIntersectingPatches")
 513                 );
 514                 patchesToRemove.transfer(patchesToRemoveCopy);
 515             }
 516
 517             const triSurf& ts = octree_.surface();
 518             boolList removeFacets(ts.size(), false);
 519
 520             //- remove facets in patches
 521             forAllConstIter(HashSet<word>, patchesToRemove, it)
 522             {
 523                 const labelList matchedPatches = ts.findPatches(it.key());
 524                 boolList activePatch(ts.patches().size(), false);
 525                 forAll(matchedPatches, ptchI)
 526                     activePatch[matchedPatches[ptchI]] = true;
 527
 528                 forAll(ts, triI)
 529                 {
 530                     if( activePatch[ts[triI].region()] )
 531                         removeFacets[triI] = true;
 532                 }
 533             }
 534
 535             //- remove facets in subsets
 536             forAllConstIter(HashSet<word>, patchesToRemove, it)
 537             {
 538                 const label subsetID = ts.facetSubsetIndex(it.key());
 539                 if( subsetID >= 0 )
 540                 {
 541                     labelLongList facets;
 542                     ts.facetsInSubset(subsetID, facets);
 543
 544                     forAll(facets, i)
 545                         removeFacets[facets[i]] = true;
 546                 }
 547             }
 548
 549             //- set BOUNDARY flag to boxes intersected by the given facets
 550             DynList<label> containedTriangles;
 551             forAll(boxType, leafI)
 552             {
 553                 octree_.containedTriangles(leafI, containedTriangles);
 554
 555                 forAll(containedTriangles, i)
 556                 {
 557                     if( removeFacets[containedTriangles[i]] )
 558                     {
 559                         boxType[leafI] = NONE;
 560                     }
 561                 }
 562             }
 563         }
 564     }
 565     else if( meshDict_.found("keepCellsIntersectingPatches") )
 566     {
 567         wordHashSet patchesToKeep;
 568
 569         if( meshDict_.isDict("keepCellsIntersectingPatches") )
 570         {
 571             const dictionary& dict =
 572                 meshDict_.subDict("keepCellsIntersectingPatches");
 573             const wordList patchNames = dict.toc();
 574
 575             forAll(patchNames, patchI)
 576                 patchesToKeep.insert(patchNames[patchI]);
 577         }
 578         else
 579         {
 580             wordHashSet patchesToKeepCopy
 581             (
 582                 meshDict_.lookup("keepCellsIntersectingPatches")
 583             );
 584             patchesToKeep.transfer(patchesToKeepCopy);
 585         }
 586
 587         const triSurf& ts = octree_.surface();
 588         boolList keepFacets(ts.size(), false);
 589
 590         //- keep facets in patches
 591         forAllConstIter(HashSet<word>, patchesToKeep, it)
 592         {
 593             const labelList matchedPatches = ts.findPatches(it.key());
 594             boolList activePatch(ts.patches().size(), false);
 595             forAll(matchedPatches, ptchI)
 596                 activePatch[matchedPatches[ptchI]] = true;
 597
 598             forAll(ts, triI)
 599             {
 600                 if( activePatch[ts[triI].region()] )
 601                     keepFacets[triI] = true;
 602             }
 603         }
 604
 605         //- keep facets in subsets
 606         forAllConstIter(wordHashSet, patchesToKeep, it)
 607         {
 608             const label subsetID = ts.facetSubsetIndex(it.key());
 609
 610             if( subsetID >= 0 )
 611             {
 612                 labelLongList facets;
 613                 ts.facetsInSubset(subsetID, facets);
 614
 615                 forAll(facets, i)
 616                     keepFacets[facets[i]] = true;
 617             }
 618         }
 619
 620         //- set MESHCELL flag to boxes intersected by the given facets
 621         DynList<label> containedTriangles;
 622         forAll(boxType, leafI)
 623         {
 624             octree_.containedTriangles(leafI, containedTriangles);
 625
 626             forAll(containedTriangles, i)
 627             {
 628                 if( keepFacets[containedTriangles[i]] )
 629                 {
 630                     boxType[leafI] = MESHCELL;
 631                 }
 632             }
 633         }
 634     }
 635
 636     //- set BOUNDARY flag to boxes which do not have a MESHCELL flag
 637     DynList<label> neighs;
 638     # ifdef USE_OMP
 639     # pragma omp parallel for if( boxType.size() > 1000 ) \
 640     private(neighs) schedule(dynamic, 20)
 641     # endif
 642     forAll(boxType, leafI)
 643     {
 644         if( boxType[leafI] & MESHCELL )
 645         {
 646             for(label i=0;i<6;++i)
 647             {
 648                 neighs.clear();
 649                 octree_.findNeighboursInDirection(leafI, i, neighs);
 650
 651                 forAll(neighs, neiI)
 652                 {
 653                     const label neiLabel = neighs[neiI];
 654
 655                     if( neiLabel < 0 )
 656                         continue;
 657
 658                     if( !(boxType[neiLabel] & MESHCELL) )
 659                         boxType[neiLabel] = BOUNDARY;
 660                 }
 661             }
 662         }
 663     }
 664
 665     if( Pstream::parRun() )
 666     {
 667         //- make sure that all processors have the same information
 668         //- about BOUNDARY boxes
 669         const labelLongList& globalLeafLabel = this->globalLeafLabel();
 670         const VRWGraph& leafAtProcs = this->leafAtProcs();
 671         const Map<label>& globalLeafToLocal =
 672             this->globalToLocalLeafAddressing();
 673
 674         std::map<label, labelLongList> exchangeData;
 675         forAll(octree_.neiProcs(), procI)
 676             exchangeData.insert
 677             (
 678                 std::make_pair
 679                 (
 680                     octree_.neiProcs()[procI],
 681                     labelLongList()
 682                 )
 683             );
 684
 685         forAllConstIter(Map<label>, globalLeafToLocal, iter)
 686         {
 687             const label leafI = iter();
 688
 689             if( boxType[leafI] & BOUNDARY )
 690             {
 691                 forAllRow(leafAtProcs, leafI, procI)
 692                 {
 693                     const label neiProc = leafAtProcs(leafI, procI);
 694
 695                     if( neiProc == Pstream::myProcNo() )
 696                         continue;
 697
 698                     exchangeData[neiProc].append(globalLeafLabel[leafI]);
 699                 }
 700             }
 701         }
 702
 703         labelLongList receivedData;
 704         help::exchangeMap(exchangeData, receivedData);
 705
 706         forAll(receivedData, i)
 707             boxType[globalLeafToLocal[receivedData[i]]] = BOUNDARY;
 708     }
 709 }
 710
 711 void meshOctreeAddressing::calculateNodeType() const
 712 {
 713     const FRWGraph<label, 8>& nodeLeaves = this->nodeLeaves();
 714
 715     nodeTypePtr_ = new List<direction>(nNodes_, NONE);
 716     List<direction>& nodeType = *nodeTypePtr_;
 717
 718     # ifdef USE_OMP
 719     # pragma omp parallel for schedule(static, 1)
 720     # endif
 721     forAll(nodeLeaves, nodeI)
 722     {
 723         forAllRow(nodeLeaves, nodeI, nlI)
 724         {
 725             const label leafI = nodeLeaves(nodeI, nlI);
 726
 727             if( leafI == -1 )
 728                 continue;
 729
 730             const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
 731
 732             if(
 733                 (oc.cubeType() & meshOctreeCubeBasic::OUTSIDE) ||
 734                 (oc.cubeType() & meshOctreeCubeBasic::UNKNOWN)
 735             )
 736             {
 737                 nodeType[nodeI] |= OUTERNODE;
 738                 break;
 739             }
 740             else if(
 741                 !octree_.hasContainedTriangles(leafI) &&
 742                 (oc.cubeType() &  meshOctreeCubeBasic::INSIDE)
 743             )
 744             {
 745                 nodeType[nodeI] |= INNERNODE;
 746                 break;
 747             }
 748         }
 749     }
 750 }
 751
 752 void meshOctreeAddressing::createOctreeFaces() const
 753 {
 754     octreeFacesPtr_ = new VRWGraph();
 755     octreeFacesOwnersPtr_ = new labelLongList();
 756     octreeFacesNeighboursPtr_ = new labelLongList();
 757
 758     const VRWGraph& nodeLabels = this->nodeLabels();
 759     const List<direction>& boxType = this->boxType();
 760     this->nodeLeaves();
 761
 762     label nFaces(0);
 763     labelList rowSizes, chunkSizes;
 764
 765     # ifdef USE_OMP
 766     # pragma omp parallel
 767     # endif
 768     {
 769         //- faces are created and stored into helper arrays, and each thread
 770         //- allocates its own graph for storing faces. The faces are generated
 771         //- by dividing the octree leaves into chunks, and distributing these
 772         //- chunks over the threads. There are four chunks per each thread to
 773         //- improve load balancing. The number of faces generated in each chunk
 774         //- is stored and later in used to store the faces into the octree faces
 775         //- graph in the correct order
 776         VRWGraph helperFaces;
 777         labelLongList helperOwner, helperNeighbour;
 778
 779         # ifdef USE_OMP
 780         const label nThreads = omp_get_num_threads();
 781         const label threadI = omp_get_thread_num();
 782         const label nChunks = 4 * omp_get_num_threads();
 783         const label chunkSize = boxType.size() / nChunks + 1;
 784         # else
 785         const label nThreads(1);
 786         const label threadI(0);
 787         const label nChunks(1);
 788         const label chunkSize = boxType.size();
 789         # endif
 790
 791         # ifdef USE_OMP
 792         # pragma omp master
 793         # endif
 794         {
 795             chunkSizes.setSize(nChunks);
 796             chunkSizes = 0;
 797         }
 798
 799         # ifdef USE_OMP
 800         # pragma omp barrier
 801         # endif
 802
 803         for
 804         (
 805             label chunkI=threadI;
 806             chunkI<nChunks;
 807             chunkI+=nThreads
 808         )
 809         {
 810             const label start = chunkSize * chunkI;
 811             const label end = Foam::min(start+chunkSize, boxType.size());
 812
 813             const label nBefore = helperFaces.size();
 814
 815             for(label leafI=start;leafI<end;++leafI)
 816             {
 817                 const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
 818
 819                 if( boxType[leafI] & MESHCELL )
 820                 {
 821                     FixedList<label, 12> edgeCentreLabel(-1);
 822                     for(label i=0;i<12;++i)
 823                         edgeCentreLabel[i] = findEdgeCentre(leafI, i);
 824
 825                     for(label fI=0;fI<6;++fI)
 826                     {
 827                         DynList<label> neighbours;
 828                         octree_.findNeighboursInDirection
 829                         (
 830                             leafI,
 831                             fI,
 832                             neighbours
 833                         );
 834
 835                         if( neighbours.size() != 1 )
 836                             continue;
 837
 838                         const label nei = neighbours[0];
 839
 840                         //- stop if the neighbour is on other processor
 841                         if( nei == meshOctreeCubeBasic::OTHERPROC )
 842                             continue;
 843
 844                         //- create face
 845                         DynList<label, 8> f;
 846                         for(label pI=0;pI<4;++pI)
 847                         {
 848                             const label nI =
 849                                 meshOctreeCubeCoordinates::faceNodes_[fI][pI];
 850                             const label feI =
 851                                 meshOctreeCubeCoordinates::faceEdges_[fI][pI];
 852
 853                             f.append(nodeLabels(leafI, nI));
 854
 855                             if( edgeCentreLabel[feI] != -1 )
 856                                 f.append(edgeCentreLabel[feI]);
 857                         }
 858
 859                         if( nei < 0 )
 860                         {
 861                             //- face is at the boundary of the octree
 862                             helperFaces.appendList(f);
 863                             helperOwner.append(leafI);
 864                             helperNeighbour.append(-1);
 865                         }
 866                         else if( boxType[nei] & MESHCELL )
 867                         {
 868                             //- face is an internal face
 869                             if( nei > leafI )
 870                             {
 871                                 helperFaces.appendList(f);
 872                                 helperOwner.append(leafI);
 873                                 helperNeighbour.append(nei);
 874                             }
 875                             else if
 876                             (
 877                                 octree_.returnLeaf(nei).level() < oc.level()
 878                             )
 879                             {
 880                                 //- append a reversed face
 881                                 label i(1);
 882                                 for(label j=f.size()-1;j>i;--j)
 883                                 {
 884                                     const label add = f[j];
 885                                     f[j] = f[i];
 886                                     f[i] = add;
 887                                     ++i;
 888                                 }
 889
 890                                 helperFaces.appendList(f);
 891                                 helperOwner.append(nei);
 892                                 helperNeighbour.append(leafI);
 893                             }
 894                         }
 895                         else if( boxType[nei] & BOUNDARY )
 896                         {
 897                             //- face is at the boundary of the mesh cells
 898                             helperFaces.appendList(f);
 899                             helperOwner.append(leafI);
 900                             helperNeighbour.append(nei);
 901                         }
 902                     }
 903                 }
 904                 else if( boxType[leafI] & BOUNDARY )
 905                 {
 906                     for(label fI=0;fI<6;++fI)
 907                     {
 908                         DynList<label> neighbours;
 909                         octree_.findNeighboursInDirection
 910                         (
 911                             leafI,
 912                             fI,
 913                             neighbours
 914                         );
 915
 916                         if( neighbours.size() != 1 )
 917                             continue;
 918                         const label nei = neighbours[0];
 919                         if( nei < 0 )
 920                             continue;
 921                         if(
 922                             (boxType[nei] & MESHCELL) &&
 923                             (octree_.returnLeaf(nei).level() < oc.level())
 924                         )
 925                         {
 926                             //- add a boundary face
 927                             const label* fNodes =
 928                                 meshOctreeCubeCoordinates::faceNodes_[fI];
 929                             face cf(4);
 930                             for(label i=0;i<4;++i)
 931                             {
 932                                 cf[i] = nodeLabels(leafI, fNodes[i]);
 933                             }
 934
 935                             helperFaces.appendList(cf.reverseFace());
 936                             helperOwner.append(nei);
 937                             helperNeighbour.append(leafI);
 938                         }
 939                     }
 940                 }
 941             }
 942
 943             //- store the size of this chunk
 944             chunkSizes[chunkI] = helperFaces.size() - nBefore;
 945         }
 946
 947         //- set the sizes of faces graph
 948         # ifdef USE_OMP
 949         # pragma omp critical
 950         # endif
 951         nFaces += helperFaces.size();
 952
 953         # ifdef USE_OMP
 954         # pragma omp barrier
 955
 956         # pragma omp master
 957         # endif
 958         {
 959             rowSizes.setSize(nFaces);
 960             octreeFacesPtr_->setSize(nFaces);
 961             octreeFacesOwnersPtr_->setSize(nFaces);
 962             octreeFacesNeighboursPtr_->setSize(nFaces);
 963         }
 964
 965         # ifdef USE_OMP
 966         # pragma omp barrier
 967         # endif
 968
 969         //- set the size of face graph rows and copy owners and neighbours
 970         for
 971         (
 972             label chunkI=threadI;
 973             chunkI<nChunks;
 974             chunkI+=nThreads
 975         )
 976         {
 977             label start(0), localStart(0);
 978             for(label i=0;i<chunkI;++i)
 979                 start += chunkSizes[i];
 980             for(label i=threadI;i<chunkI;i+=nThreads)
 981                 localStart += chunkSizes[i];
 982
 983             for(label faceI=0;faceI<chunkSizes[chunkI];++faceI)
 984             {
 985                 octreeFacesOwnersPtr_->operator[](start) =
 986                     helperOwner[localStart];
 987                 octreeFacesNeighboursPtr_->operator[](start) =
 988                     helperNeighbour[localStart];
 989                 rowSizes[start++] = helperFaces.sizeOfRow(localStart++);
 990             }
 991         }
 992
 993         # ifdef USE_OMP
 994         # pragma omp barrier
 995
 996         //- set the size of octree faces
 997         # pragma omp master
 998         # endif
 999         VRWGraphSMPModifier(*octreeFacesPtr_).setSizeAndRowSize(rowSizes);
1000
1001         # ifdef USE_OMP
1002         # pragma omp barrier
1003         # endif
1004
1005         //- copy the data into octree faces
1006         for
1007         (
1008             label chunkI=threadI;
1009             chunkI<nChunks;
1010             chunkI+=nThreads
1011         )
1012         {
1013             label start(0), localStart(0);
1014
1015             for(label i=0;i<chunkI;++i)
1016                 start += chunkSizes[i];
1017             for(label i=threadI;i<chunkI;i+=nThreads)
1018                 localStart += chunkSizes[i];
1019
1020             for(label faceI=0;faceI<chunkSizes[chunkI];++faceI)
1021             {
1022                 for(label i=0;i<helperFaces.sizeOfRow(localStart);++i)
1023                     octreeFacesPtr_->operator()(start, i) =
1024                         helperFaces(localStart, i);
1025
1026                 ++start;
1027                 ++localStart;
1028             }
1029         }
1030     }
1031
1032     # ifdef DEBUGVrt
1033     List<vector> sum(octree_.numberOfLeaves(), vector::zero);
1034     for(label faceI=0;faceI<octreeFacesPtr_->size();++faceI)
1035     {
1036         face f(octreeFacesPtr_->sizeOfRow(faceI));
1037         forAll(f, pI)
1038             f[pI] = octreeFacesPtr_->operator()(faceI, pI);
1039         const vector n = f.normal(this->octreePoints());
1040
1041         sum[(*octreeFacesOwnersPtr_)[faceI]] += n;
1042         const label nei = (*octreeFacesNeighboursPtr_)[faceI];
1043
1044         if( nei < 0 )
1045             continue;
1046         sum[nei] -= n;
1047     }
1048
1049     forAll(sum, lfI)
1050     {
1051         if
1052         (
1053             Pstream::parRun() &&
1054             octree_.returnLeaf(lfI).procNo() != Pstream::myProcNo()
1055         )
1056             continue;
1057         if( (boxType[lfI] & MESHCELL) && (mag(sum[lfI]) > SMALL) )
1058             Info << "Leaf " << lfI << " is not closed " << sum[lfI] << endl;
1059     }
1060     # endif
1061 }
1062
1063 void meshOctreeAddressing::calculateLeafFaces() const
1064 {
1065     const labelLongList& owner = octreeFaceOwner();
1066     const labelLongList& neighbour = octreeFaceNeighbour();
1067
1068     leafFacesPtr_ = new VRWGraph(octree_.numberOfLeaves());
1069     VRWGraph& leafFaces = *leafFacesPtr_;
1070
1071     labelList nlf(leafFaces.size(), 0);
1072     forAll(owner, fI)
1073     {
1074         ++nlf[owner[fI]];
1075         if( neighbour[fI] < 0 )
1076             continue;
1077         ++nlf[neighbour[fI]];
1078     }
1079
1080     forAll(nlf, leafI)
1081         leafFaces.setRowSize(leafI, nlf[leafI]);
1082     nlf = 0;
1083
1084     forAll(owner, fI)
1085     {
1086         leafFaces(owner[fI], nlf[owner[fI]]++) = fI;
1087         if( neighbour[fI] < 0 )
1088             continue;
1089         leafFaces(neighbour[fI], nlf[neighbour[fI]]++) = fI;
1090     }
1091 }
1092
1093 void meshOctreeAddressing::calculateNodeFaces() const
1094 {
1095     const VRWGraph& octreeFaces = this->octreeFaces();
1096     nodeFacesPtr_ = new VRWGraph(numberOfNodes());
1097     VRWGraph& nodeFaces = *nodeFacesPtr_;
1098
1099     VRWGraphSMPModifier(nodeFaces).reverseAddressing(octreeFaces);
1100     nodeFaces.setSize(numberOfNodes());
1101 }
1102
1103 void meshOctreeAddressing::calculateLeafLeaves() const
1104 {
1105     const labelLongList& owner = octreeFaceOwner();
1106     const labelLongList& neighbour = octreeFaceNeighbour();
1107
1108     leafLeavesPtr_ = new VRWGraph(octree_.numberOfLeaves());
1109     VRWGraph& leafLeaves = *leafLeavesPtr_;
1110
1111     labelList nNei(leafLeaves.size(), 0);
1112     forAll(owner, faceI)
1113     {
1114         if( owner[faceI] < 0 )
1115             continue;
1116         if( neighbour[faceI] < 0 )
1117             continue;
1118
1119         ++nNei[owner[faceI]];
1120         ++nNei[neighbour[faceI]];
1121     }
1122
1123     forAll(nNei, leafI)
1124         leafLeaves.setRowSize(leafI, nNei[leafI]);
1125
1126     nNei = 0;
1127
1128     forAll(owner, faceI)
1129     {
1130         if( owner[faceI] < 0 )
1131             continue;
1132         if( neighbour[faceI] < 0 )
1133             continue;
1134
1135         leafLeaves(owner[faceI], nNei[owner[faceI]]++) = neighbour[faceI];
1136         leafLeaves(neighbour[faceI], nNei[neighbour[faceI]]++) = owner[faceI];
1137     }
1138 }
1139
1140 void meshOctreeAddressing::createOctreeEdges() const
1141 {
1142     const VRWGraph& faces = this->octreeFaces();
1143
1144     //- allocate memory for edges, face-edges addressing
1145     //- and node-edges addressing
1146     octreeEdgesPtr_ = new LongList<edge>();
1147     LongList<edge>& edges = *octreeEdgesPtr_;
1148     faceEdgesPtr_ = new VRWGraph(faces.size());
1149     VRWGraph& faceEdges = *faceEdgesPtr_;
1150     nodeEdgesPtr_ = new VRWGraph();
1151     VRWGraph& nodeEdges = *nodeEdgesPtr_;
1152     nodeEdges.setSizeAndColumnWidth(nNodes_, 6);
1153
1154     forAll(faces, faceI)
1155     {
1156         faceEdges.setRowSize(faceI, faces[faceI].size());
1157         forAllRow(faceEdges, faceI, feI)
1158             faceEdges(faceI, feI) = -1;
1159     }
1160
1161     forAll(faces, faceI)
1162     {
1163         const label nEdges = faces.sizeOfRow(faceI);
1164
1165         for(label eI=0;eI<nEdges;++eI)
1166         {
1167             const edge e
1168             (
1169                 faces(faceI, eI),
1170                 faces(faceI, (eI+1)%nEdges)
1171             );
1172
1173             label eLabel(-1);
1174             forAllRow(nodeEdges, e.start(), neI)
1175             {
1176                 if( edges[nodeEdges(e.start(), neI)] == e )
1177                 {
1178                     eLabel = nodeEdges(e.start(), neI);
1179                     break;
1180                 }
1181             }
1182
1183             if( eLabel < 0 )
1184             {
1185                 //- append new edge
1186                 faceEdges(faceI, eI) = edges.size();
1187                 nodeEdges.append(e.start(), edges.size());
1188                 nodeEdges.append(e.end(), edges.size());
1189
1190                 edges.append(e);
1191             }
1192             else
1193             {
1194                 faceEdges(faceI, eI) = eLabel;
1195             }
1196         }
1197     }
1198 }
1199
1200 void meshOctreeAddressing::calculateLeafEdges() const
1201 {
1202     const VRWGraph& edgeLeaves = this->edgeLeaves();
1203
1204     leafEdgesPtr_ = new VRWGraph();
1205     VRWGraph& leafEdges = *leafEdgesPtr_;
1206
1207     VRWGraphSMPModifier(leafEdges).reverseAddressing(edgeLeaves);
1208     leafEdges.setSize(octree_.numberOfLeaves());
1209 }
1210
1211 void meshOctreeAddressing::calculateEdgeLeaves() const
1212 {
1213     const VRWGraph& edgeFaces = this->edgeFaces();
1214     const labelLongList& owner = this->octreeFaceOwner();
1215     const labelLongList& neighbour = this->octreeFaceNeighbour();
1216
1217     edgeLeavesPtr_ = new VRWGraph();
1218     VRWGraph& edgeLeaves = *edgeLeavesPtr_;
1219     edgeLeaves.setSizeAndColumnWidth(edgeFaces.size(), 4);
1220
1221     forAll(edgeFaces, edgeI)
1222     {
1223         forAllRow(edgeFaces, edgeI, efI)
1224         {
1225             const label fI = edgeFaces(edgeI, efI);
1226             const label own = owner[fI];
1227             const label nei = neighbour[fI];
1228
1229             edgeLeaves.appendIfNotIn(edgeI, own);
1230
1231             if( nei < 0 )
1232                 continue;
1233             edgeLeaves.appendIfNotIn(edgeI, nei);
1234         }
1235     }
1236 }
1237
1238 void meshOctreeAddressing::calculateEdgeFaces() const
1239 {
1240     const VRWGraph& faceEdges = this->faceEdges();
1241     edgeFacesPtr_ = new VRWGraph(octreeEdges().size());
1242     VRWGraph& edgeFaces = *edgeFacesPtr_;
1243
1244     VRWGraphSMPModifier(edgeFaces).reverseAddressing(faceEdges);
1245     edgeFaces.setSize(octreeEdges().size());
1246 }
1247
1248 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1249
1250 } // End namespace Foam
1251
1252 // ************************************************************************* //