src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctree.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 "meshOctree.H"
  29 #include "triSurf.H"
  30 #include "boundBox.H"
  31 #include "demandDrivenData.H"
  32
  33 # ifdef USE_OMP
  34 #include <omp.h>
  35 # endif
  36
  37 //#define DEBUGSearch
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 // Construct from surface
  47 meshOctree::meshOctree(const triSurf& ts, const bool isQuadtree)
  48 :
  49     surface_(ts),
  50     neiProcs_(),
  51     neiRange_(),
  52     initialCubePtr_(NULL),
  53     initialCubeRotation_(0),
  54     rootBox_(),
  55     isRootInitialised_(false),
  56     searchRange_(0.0),
  57     octantVectors_(),
  58     vrtLeavesPos_(),
  59     regularityPositions_(),
  60     dataSlots_(),
  61     leaves_(),
  62     isQuadtree_(isQuadtree)
  63 {
  64     createInitialOctreeBox();
  65
  66     setOctantVectorsAndPositions();
  67 }
  68
  69 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  70
  71 meshOctree::~meshOctree()
  72 {}
  73
  74 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  75
  76 void meshOctree::setOctantVectorsAndPositions()
  77 {
  78     octantVectors_[0] = Vector<label>(-1, -1, -1);
  79     octantVectors_[1] = Vector<label>(1, -1, -1);
  80     octantVectors_[2] = Vector<label>(-1, 1, -1);
  81     octantVectors_[3] = Vector<label>(1, 1, -1);
  82     octantVectors_[4] = Vector<label>(-1, -1, 1);
  83     octantVectors_[5] = Vector<label>(1, -1, 1);
  84     octantVectors_[6] = Vector<label>(-1, 1, 1);
  85     octantVectors_[7] = Vector<label>(1, 1, 1);
  86
  87     # ifdef DEBUGSearch
  88     Info << "Octant vectors " << octantVectors_ << endl;
  89     # endif
  90
  91     //- set regularity positions
  92     //- neighbours over faces
  93     regularityPositions_[0] = meshOctreeCubeCoordinates(-1, 0, 0, 0);
  94     regularityPositions_[1] = meshOctreeCubeCoordinates(1, 0, 0, 0);
  95     regularityPositions_[2] = meshOctreeCubeCoordinates(0, -1, 0, 0);
  96     regularityPositions_[3] = meshOctreeCubeCoordinates(0, 1, 0, 0);
  97     regularityPositions_[4] = meshOctreeCubeCoordinates(0, 0, -1, 0);
  98     regularityPositions_[5] = meshOctreeCubeCoordinates(0, 0, 1, 0);
  99
 100     //- neighbours over edges
 101     //- edges in x-direction
 102     regularityPositions_[6] = meshOctreeCubeCoordinates(0, -1, -1, 0);
 103     regularityPositions_[7] = meshOctreeCubeCoordinates(0, 1, -1, 0);
 104     regularityPositions_[8] = meshOctreeCubeCoordinates(0, -1, 1, 0);
 105     regularityPositions_[9] = meshOctreeCubeCoordinates(0, 1, 1, 0);
 106
 107     //- edges in y-direction
 108     regularityPositions_[10] = meshOctreeCubeCoordinates(-1, 0, -1, 0);
 109     regularityPositions_[11] = meshOctreeCubeCoordinates(1, 0, -1, 0);
 110     regularityPositions_[12] = meshOctreeCubeCoordinates(-1, 0, 1, 0);
 111     regularityPositions_[13] = meshOctreeCubeCoordinates(1, 0, 1, 0);
 112
 113     //- edges in z-direction
 114     regularityPositions_[14] = meshOctreeCubeCoordinates(-1, -1, 0, 0);
 115     regularityPositions_[15] = meshOctreeCubeCoordinates(1, -1, 0, 0);
 116     regularityPositions_[16] = meshOctreeCubeCoordinates(-1, 1, 0, 0);
 117     regularityPositions_[17] = meshOctreeCubeCoordinates(1, 1, 0, 0);
 118
 119     //- neighbours over vertices
 120     regularityPositions_[18] = meshOctreeCubeCoordinates(-1, -1, -1, 0);
 121     regularityPositions_[19] = meshOctreeCubeCoordinates(1, -1, -1, 0);
 122     regularityPositions_[20] = meshOctreeCubeCoordinates(-1, 1, -1, 0);
 123     regularityPositions_[21] = meshOctreeCubeCoordinates(1, 1, -1, 0);
 124     regularityPositions_[22] = meshOctreeCubeCoordinates(-1, -1, 1, 0);
 125     regularityPositions_[23] = meshOctreeCubeCoordinates(1, -1, 1, 0);
 126     regularityPositions_[24] = meshOctreeCubeCoordinates(-1, 1, 1, 0);
 127     regularityPositions_[25] = meshOctreeCubeCoordinates(1, 1, 1, 0);
 128
 129     # ifdef DEBUGSearch
 130     Info << "Regularity positions " << regularityPositions_ << endl;
 131     # endif
 132
 133     //- set vrtLeavesPos_
 134     for(label vrtI=0;vrtI<8;++vrtI)
 135     {
 136         FixedList<label, 3> vc(0);
 137
 138         if( vrtI & 1 )
 139             vc[0] += 1;
 140         if( vrtI & 2 )
 141             vc[1] += 1;
 142         if( vrtI & 4 )
 143             vc[2] += 1;
 144
 145         # ifdef DEBUGSearch
 146         Info << "Vert " << vrtI << " vc " << vc << endl;
 147         # endif
 148
 149         for(label i=0;i<8;++i)
 150         {
 151             FixedList<label, 3> pos;
 152
 153             for(label j=0;j<3;++j)
 154             {
 155                 if( vc[j] == 0 && octantVectors_[i][j] == 1 )
 156                 {
 157                     pos[j] = 0;
 158                 }
 159                 else if( vc[j] == 1 && octantVectors_[i][j] == 1 )
 160                 {
 161                     pos[j] = 1;
 162                 }
 163                 else
 164                 {
 165                     pos[j] = vc[j] + octantVectors_[i][j];
 166                 }
 167             }
 168
 169             vrtLeavesPos_[vrtI][i] =
 170                 meshOctreeCubeCoordinates
 171                 (
 172                     pos[0],
 173                     pos[1],
 174                     pos[2],
 175                     direction(0)
 176                 );
 177         }
 178     }
 179
 180     # ifdef DEBUGSearch
 181     Info << "vrtLeavesPos_ " << vrtLeavesPos_ << endl;
 182     # endif
 183 }
 184
 185 void meshOctree::createInitialOctreeBox()
 186 {
 187     //- create initial octree box
 188     boundBox bb(surface_.points());
 189     const point& min_ = bb.min();
 190     const point& max_ = bb.max();
 191
 192     const point c = (max_ + min_) / 2.0;
 193     scalar cs = 1.5 * (max_.x() - min_.x()) / 2.0;
 194     if( cs < (1.5 * (max_.y() - min_.y()) / 2.0) )
 195     {
 196         cs = 1.5 * (max_.y() - min_.y()) / 2.0;
 197     }
 198     if( cs < (1.5 * (max_.z() - min_.z()) / 2.0) )
 199     {
 200         cs = 1.5 * (max_.z() - min_.z()) / 2.0;
 201     }
 202
 203     //- create root box and initial cube
 204     rootBox_ = boundBox(c - point(cs, cs, cs), c + point(cs, cs, cs));
 205
 206     if( Pstream::parRun() )
 207     {
 208         reduce(rootBox_.min(), minOp<point>());
 209         reduce(rootBox_.max(), maxOp<point>());
 210     }
 211
 212     //- allocate data slots
 213     # ifdef USE_OMP
 214     if( omp_get_num_procs() > 0 )
 215     {
 216         dataSlots_.setSize(omp_get_num_procs());
 217     }
 218     else
 219     {
 220         dataSlots_.setSize(1);
 221     }
 222     # else
 223     dataSlots_.setSize(1);
 224     # endif
 225
 226     meshOctreeSlot* slotPtr = &dataSlots_[0];
 227
 228     if( !isQuadtree_ )
 229     {
 230         slotPtr->cubes_.append
 231         (
 232             meshOctreeCube
 233             (
 234                 meshOctreeCubeCoordinates(0, 0, 0, 0),
 235                 surface_.size(),
 236                 slotPtr
 237             )
 238         );
 239     }
 240     else
 241     {
 242         slotPtr->cubes_.append
 243         (
 244             meshOctreeCube
 245             (
 246                 meshOctreeCubeCoordinates(0, 0, -10, 0),
 247                 surface_.size(),
 248                 slotPtr
 249             )
 250         );
 251     }
 252
 253     initialCubePtr_ = &slotPtr->cubes_[0];
 254
 255     leaves_.setSize(1);
 256     leaves_[0] = initialCubePtr_;
 257 }
 258
 259 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 260
 261 } // End namespace Foam
 262
 263 // ************************************************************************* //