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