src/meshLibrary/utilities/octrees/meshOctree/meshOctreeModifier/meshOctreeModifierRefineSelectedBoxes.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 "meshOctreeModifier.H"
  29 #include "triSurf.H"
  30 #include "HashSet.H"
  31
  32 # ifdef USE_OMP
  33 #include <omp.h>
  34 # endif
  35
  36 #include <sys/stat.h>
  37
  38 //#define OCTREETiming
  39 //#define DEBUGSearch
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 void meshOctreeModifier::markAdditionalLayers
  49 (
  50     List<direction>& refineBox,
  51     const direction nLayers
  52 ) const
  53 {
  54     const FixedList<meshOctreeCubeCoordinates, 26>& rp =
  55         octree_.regularityPositions_;
  56     const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
  57
  58     //- this is needed for parallel runs to reduce the bandwidth
  59     labelHashSet transferCoordinates;
  60
  61     FixedList<meshOctreeCube*, 26> neighbours;
  62
  63     for(direction i=1;i<=nLayers;++i)
  64     {
  65         LongList<meshOctreeCubeCoordinates> processorChecks;
  66
  67         # ifdef USE_OMP
  68         # pragma omp parallel for if( leaves.size() > 1000 ) \
  69         private(neighbours) schedule(dynamic, 20)
  70         # endif
  71         forAll(leaves, leafI)
  72         {
  73             if( refineBox[leafI] != i )
  74                 continue;
  75
  76             const meshOctreeCube* oc = leaves[leafI];
  77
  78             forAll(rp, posI)
  79             {
  80                 const meshOctreeCubeCoordinates cc
  81                 (
  82                     oc->coordinates() + rp[posI]
  83                 );
  84
  85                 const label neiLabel = octree_.findLeafLabelForPosition(cc);
  86
  87                 if( neiLabel > -1 )
  88                 {
  89                     neighbours[posI] = leaves[neiLabel];
  90                 }
  91                 else if( neiLabel == -1 )
  92                 {
  93                     neighbours[posI] = NULL;
  94                 }
  95                 else if( neiLabel == meshOctreeCubeBasic::OTHERPROC )
  96                 {
  97                     neighbours[posI] = NULL;
  98
  99                     # ifdef USE_OMP
 100                     # pragma omp critical
 101                     # endif
 102                     {
 103                         if( !transferCoordinates.found(leafI) )
 104                         {
 105                             processorChecks.append(oc->coordinates());
 106                             transferCoordinates.insert(leafI);
 107                         }
 108                     }
 109                 }
 110             }
 111
 112             forAll(neighbours, neiI)
 113             {
 114                 const meshOctreeCube* nei = neighbours[neiI];
 115                 if( !nei ) continue;
 116                 if( !nei->isLeaf() ) continue;
 117                 if( nei->level() > oc->level() ) continue;
 118
 119                 if( !refineBox[nei->cubeLabel()] )
 120                 {
 121                     refineBox[nei->cubeLabel()] = i+1;
 122                 }
 123             }
 124         }
 125
 126         if( octree_.neiProcs().size() )
 127         {
 128             LongList<meshOctreeCubeCoordinates> receivedCoords;
 129             octree_.exchangeRequestsWithNeighbourProcessors
 130             (
 131                 processorChecks,
 132                 receivedCoords
 133             );
 134
 135             //- check consistency with received cube coordinates
 136             # ifdef USE_OMP
 137             # pragma omp parallel for if( receivedCoords.size() > 1000 ) \
 138             schedule(dynamic, 20)
 139             # endif
 140             forAll(receivedCoords, ccI)
 141             {
 142                 forAll(rp, posI)
 143                 {
 144                     const meshOctreeCubeCoordinates cc
 145                     (
 146                         receivedCoords[ccI] + rp[posI]
 147                     );
 148
 149                     const meshOctreeCube* nei =
 150                         octree_.findCubeForPosition(cc);
 151
 152                     if( !nei ) continue;
 153                     if( !nei->isLeaf() ) continue;
 154                     if( nei->level() > cc.level() ) continue;
 155
 156                     if( !refineBox[nei->cubeLabel()] )
 157                     {
 158                         refineBox[nei->cubeLabel()] = i+1;
 159                     }
 160                 }
 161             }
 162         }
 163     }
 164 }
 165
 166 void meshOctreeModifier::refineSelectedBoxes
 167 (
 168     List<direction>& refineBox,
 169     const bool hexRefinement
 170 )
 171 {
 172     # ifdef OCTREETiming
 173     const scalar startTime = omp_get_wtime();
 174     # endif
 175
 176     //- ensure that refinement will produce 1-irregular octree
 177     do
 178     {
 179         ensureCorrectRegularity(refineBox);
 180     } while( hexRefinement && ensureCorrectRegularitySons(refineBox) );
 181
 182     # ifdef OCTREETiming
 183     const scalar regTime = omp_get_wtime();
 184     Info << "Time for ensuring regularity " << (regTime-startTime) << endl;
 185     # endif
 186
 187     const triSurf& surface = octree_.surface();
 188     const boundBox& rootBox = octree_.rootBox();
 189     const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
 190
 191     //- this is needed for thread safety
 192     //- such solutions make me a sad bunny :(
 193     surface.facetEdges();
 194     surface.edgeFacets();
 195     surface.edges();
 196
 197     # ifdef USE_OMP
 198     # pragma omp parallel num_threads(octree_.dataSlots_.size())
 199     # endif
 200     {
 201         # ifdef USE_OMP
 202         meshOctreeSlot* slotPtr = &octree_.dataSlots_[omp_get_thread_num()];
 203         # else
 204         meshOctreeSlot* slotPtr = &octree_.dataSlots_[0];
 205         # endif
 206
 207         if( !octree_.isQuadtree() )
 208         {
 209             # ifdef USE_OMP
 210             # pragma omp for schedule(dynamic, 100)
 211             # endif
 212             forAll(leaves, leafI)
 213             {
 214                 if( refineBox[leafI] )
 215                     leaves[leafI]->refineCube(surface, rootBox, slotPtr);
 216             }
 217         }
 218         else
 219         {
 220             # ifdef USE_OMP
 221             # pragma omp for schedule(dynamic, 100)
 222             # endif
 223             forAll(leaves, leafI)
 224             {
 225                 if( refineBox[leafI] )
 226                     leaves[leafI]->refineCube2D(surface, rootBox, slotPtr);
 227             }
 228         }
 229     }
 230
 231     createListOfLeaves();
 232
 233     # ifdef OCTREETiming
 234     Info << "Time for actual refinement " << (omp_get_wtime()-regTime) << endl;
 235     # endif
 236 }
 237
 238 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 239
 240 } // End namespace Foam
 241
 242 // ************************************************************************* //