src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeNeighbourSearches.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 "boundBox.H"
  30
  31 //#define OCTREE_DEBUG
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39 // Private member functions
  40
  41 label meshOctree::findLeafContainingVertex
  42 (
  43     const point& p
  44 ) const
  45 {
  46     # ifdef OCTREE_DEBUG
  47     Info << "Finding leaf for vertex " << p << endl;
  48     # endif
  49
  50     const meshOctreeCube* ocPtr = initialCubePtr_;
  51
  52     if( !ocPtr->isVertexInside(rootBox_, p) )
  53     {
  54         # ifdef OCTREE_DEBUG
  55         Info << "Vertex " << p << " is not in the initial cube" << endl;
  56         # endif
  57
  58         FatalErrorIn
  59         (
  60             "label meshOctree::findLeafContainingVertex(const point&) const"
  61         ) << "Point " << p << " is not inside the initial cube" << endl;
  62
  63         throw "Found invalid locations of points";
  64
  65         return -1;
  66     }
  67
  68     bool finished(false);
  69
  70     do
  71     {
  72         if( ocPtr && !ocPtr->isLeaf() )
  73         {
  74             //- find a subCube containing the vertex;
  75             const point c = ocPtr->centre(rootBox_);
  76
  77             label scI(0);
  78
  79             if( p.x() >= c.x() )
  80                 scI |= 1;
  81             if( p.y() >= c.y() )
  82                 scI |= 2;
  83             if( !isQuadtree_ && p.z() >= c.z() )
  84                 scI |= 4;
  85
  86             ocPtr = ocPtr->subCube(scI);
  87         }
  88         else
  89         {
  90             finished = true;
  91         }
  92     } while( !finished );
  93
  94     if( ocPtr )
  95     {
  96         return ocPtr->cubeLabel();
  97     }
  98
  99     return meshOctreeCubeBasic::OTHERPROC;
 100 }
 101
 102 void meshOctree::findLeavesInSphere
 103 (
 104     const point& c,
 105     const scalar r,
 106     DynList<label>& containedLeaves
 107 ) const
 108 {
 109     containedLeaves.clear();
 110
 111     initialCubePtr_->leavesInSphere(rootBox_, c, r, containedLeaves);
 112 }
 113
 114 label meshOctree::findNeighbourOverNode
 115 (
 116     const meshOctreeCubeCoordinates& cc,
 117     const label nodeI
 118 ) const
 119 {
 120     if( isQuadtree_ )
 121         return -1;
 122
 123     const meshOctreeCubeCoordinates nc(cc + regularityPositions_[18+nodeI]);
 124
 125     const meshOctreeCube* neiPtr = findCubeForPosition(nc);
 126
 127     if( !neiPtr )
 128     {
 129         const label levelLimiter = (1 << cc.level());
 130         if(
 131             (nc.posX() >= levelLimiter) || (nc.posX() < 0) ||
 132             (nc.posY() >= levelLimiter) || (nc.posY() < 0) ||
 133             (!isQuadtree_ && (nc.posZ() >= levelLimiter || nc.posZ() < 0)) ||
 134             (isQuadtree_ && (nc.posZ() != initialCubePtr_->posZ()))
 135         )
 136         {
 137             return -1;
 138         }
 139         else if( Pstream::parRun() )
 140         {
 141             return meshOctreeCubeBasic::OTHERPROC;
 142         }
 143
 144         return -1;
 145     }
 146     else if( neiPtr->isLeaf() )
 147     {
 148         # ifdef OCTREE_DEBUG
 149         if( leaves_[neiPtr->cubeLabel()] != neiPtr )
 150             FatalError << "Cube does not correspond to itself"
 151                 << abort(FatalError);
 152         # endif
 153         return neiPtr->cubeLabel();
 154     }
 155     else
 156     {
 157         FixedList<label, 8> sc(-1);
 158         for(label scI=0;scI<8;++scI)
 159         {
 160             meshOctreeCube* scPtr = neiPtr->subCube(scI);
 161             if( scPtr )
 162             {
 163                 sc[scI] = scPtr->cubeLabel();
 164             }
 165             else if( Pstream::parRun() )
 166             {
 167                 sc[scI] = meshOctreeCubeBasic::OTHERPROC;
 168             }
 169         }
 170
 171         return sc[7-nodeI];
 172     }
 173
 174     FatalErrorIn
 175     (
 176         "label meshOctree::findNeighbourOverNode("
 177         "const meshOctreeCubeCoordinates& cc,"
 178         "const label nodeI) const"
 179     ) << "Should not be here!" << abort(FatalError);
 180
 181     return -1;
 182 }
 183
 184 void meshOctree::findNeighboursOverEdge
 185 (
 186     const meshOctreeCubeCoordinates& cc,
 187     const label eI,
 188     DynList<label>& neighbourLeaves
 189 ) const
 190 {
 191     if( isQuadtree_ && (eI < 8) )
 192     {
 193         neighbourLeaves.append(-1);
 194         return;
 195     }
 196
 197     const meshOctreeCubeCoordinates nc(cc + regularityPositions_[6+eI]);
 198
 199     const meshOctreeCube* neiPtr = findCubeForPosition(nc);
 200
 201     if( !neiPtr )
 202     {
 203         const label levelLimiter = (1 << cc.level());
 204         if(
 205             (nc.posX() >= levelLimiter) || (nc.posX() < 0) ||
 206             (nc.posY() >= levelLimiter) || (nc.posY() < 0) ||
 207             (!isQuadtree_ && (nc.posZ() >= levelLimiter || nc.posZ() < 0)) ||
 208             (isQuadtree_ && (nc.posZ() != initialCubePtr_->posZ()))
 209         )
 210         {
 211             neighbourLeaves.append(-1);
 212         }
 213         else if( Pstream::parRun() )
 214         {
 215             neighbourLeaves.append(meshOctreeCubeBasic::OTHERPROC);
 216         }
 217     }
 218     else if( neiPtr->isLeaf() )
 219     {
 220         # ifdef OCTREE_DEBUG
 221         if( leaves_[neiPtr->cubeLabel()] != neiPtr )
 222             FatalError << "Cube does not correspond to itself"
 223                 << abort(FatalError);
 224         # endif
 225         neighbourLeaves.append(neiPtr->cubeLabel());
 226     }
 227     else
 228     {
 229         FixedList<label, 8> sc(-1);
 230         for(label scI=0;scI<8;++scI)
 231         {
 232             meshOctreeCube* scPtr = neiPtr->subCube(scI);
 233
 234             if( scPtr )
 235             {
 236                 sc[scI] = scPtr->cubeLabel();
 237             }
 238             else if( Pstream::parRun() )
 239             {
 240                 sc[scI] = meshOctreeCubeBasic::OTHERPROC;
 241             }
 242         }
 243
 244         const label* eNodes = meshOctreeCubeCoordinates::edgeNodes_[eI];
 245
 246         if( !isQuadtree_)
 247         {
 248             neighbourLeaves.append(sc[7-eNodes[1]]);
 249             neighbourLeaves.append(sc[7-eNodes[0]]);
 250         }
 251         else
 252         {
 253             if( sc[7-eNodes[1]] >= 0 )
 254                 neighbourLeaves.append(sc[7-eNodes[1]]);
 255             if( (sc[7-eNodes[0]] >= 0) && (sc[7-eNodes[0]] != sc[7-eNodes[1]]) )
 256                 neighbourLeaves.append(sc[7-eNodes[0]]);
 257         }
 258     }
 259 }
 260
 261 void meshOctree::findNeighboursInDirection
 262 (
 263     const meshOctreeCubeCoordinates& cc,
 264     const label dir,
 265     DynList<label>& neighbourLeaves
 266 ) const
 267 {
 268     if( isQuadtree_ && dir >= 4 )
 269     {
 270         neighbourLeaves.append(-1);
 271         return;
 272     }
 273
 274     label cpx = cc.posX();
 275     label cpy = cc.posY();
 276     label cpz = cc.posZ();
 277     switch( dir )
 278     {
 279         case 0:
 280         {
 281             cpx -= 1;
 282         }
 283         break;
 284         case 1:
 285         {
 286             cpx += 1;
 287         }
 288         break;
 289         case 2:
 290         {
 291             cpy -= 1;
 292         }
 293         break;
 294         case 3:
 295         {
 296             cpy += 1;
 297         }
 298         break;
 299         case 4:
 300         {
 301             cpz -= 1;
 302         }
 303         break;
 304         case 5:
 305         {
 306             cpz += 1;
 307         }
 308         break;
 309     }
 310
 311     const meshOctreeCube* neiPtr =
 312         findCubeForPosition
 313         (
 314             meshOctreeCubeCoordinates(cpx, cpy, cpz, cc.level())
 315         );
 316
 317     if( !neiPtr )
 318     {
 319         const label levelLimiter = (1 << cc.level());
 320         if(
 321             (cpx >= levelLimiter) || (cpx < 0) ||
 322             (cpy >= levelLimiter) || (cpy < 0) ||
 323             (!isQuadtree_ && (cpz >= levelLimiter || cpz < 0)) ||
 324             (isQuadtree_ && (cpz != initialCubePtr_->posZ()))
 325         )
 326         {
 327             neighbourLeaves.append(-1);
 328         }
 329         else if( Pstream::parRun() )
 330         {
 331             neighbourLeaves.append(meshOctreeCubeBasic::OTHERPROC);
 332         }
 333     }
 334     else if( neiPtr->isLeaf() )
 335     {
 336         # ifdef OCTREE_DEBUG
 337         if( leaves_[neiPtr->cubeLabel()] != neiPtr )
 338             FatalError << "Cube does not correspond to itself"
 339                 << abort(FatalError);
 340         # endif
 341         neighbourLeaves.append(neiPtr->cubeLabel());
 342     }
 343     else
 344     {
 345         FixedList<label, 8> sc(-1);
 346         for(label scI=0;scI<8;++scI)
 347         {
 348             meshOctreeCube* scPtr = neiPtr->subCube(scI);
 349
 350             if( scPtr )
 351             {
 352                 sc[scI] = scPtr->cubeLabel();
 353             }
 354             else if( Pstream::parRun() )
 355             {
 356                 sc[scI] = meshOctreeCubeBasic::OTHERPROC;
 357             }
 358         }
 359
 360         const label* fNodes = meshOctreeCubeCoordinates::faceNodes_[dir];
 361         for(label i=0;i<4;++i)
 362         {
 363             if( isQuadtree_ && sc[7-fNodes[i]] < 0 )
 364                 continue;
 365
 366             neighbourLeaves.append(sc[7-fNodes[i]]);
 367         }
 368     }
 369 }
 370
 371 void meshOctree::findNeighboursForLeaf
 372 (
 373     const meshOctreeCubeCoordinates& cc,
 374     DynList<label>& neighbourLeaves
 375 ) const
 376 {
 377     neighbourLeaves.clear();
 378
 379     const label nCubeFaces = isQuadtree_?4:6;
 380     for(label i=0;i<nCubeFaces;++i)
 381     {
 382         findNeighboursInDirection(cc, i, neighbourLeaves);
 383     }
 384 }
 385
 386 void meshOctree::findAllLeafNeighbours
 387 (
 388     const meshOctreeCubeCoordinates& cc,
 389     DynList<label>& neighbourLeaves
 390 ) const
 391 {
 392     neighbourLeaves.clear();
 393
 394     //- neighbours over nodes
 395     if( !isQuadtree_ )
 396     {
 397         for(label i=0;i<8;++i)
 398             neighbourLeaves.append(findNeighbourOverNode(cc, i));
 399
 400         //- neighbours over edges
 401         for(label i=0;i<12;++i)
 402             findNeighboursOverEdge(cc, i, neighbourLeaves);
 403
 404         //- neighbours over faces
 405         for(label i=0;i<6;++i)
 406             findNeighboursInDirection(cc, i, neighbourLeaves);
 407     }
 408     else
 409     {
 410         //- neighbours of an quadtree leaf
 411         //- neighbours over edges
 412         for(label i=8;i<12;++i)
 413             findNeighboursOverEdge(cc, i, neighbourLeaves);
 414
 415         //- neighbours over faces
 416         for(label i=0;i<4;++i)
 417             findNeighboursInDirection(cc, i, neighbourLeaves);
 418     }
 419 }
 420
 421 void meshOctree::findLeavesForCubeVertex
 422 (
 423     const label leafI,
 424     const direction vrtI,
 425     FixedList<label, 8>& neighbours
 426 ) const
 427 {
 428     const meshOctreeCube* oc = leaves_[leafI];
 429     const meshOctreeCubeCoordinates cc = oc->refineForPosition(vrtI);
 430
 431     FixedList<meshOctreeCubeCoordinates, 8> positions;
 432
 433     for(label i=0;i<8;++i)
 434     {
 435         positions[i] = cc + vrtLeavesPos_[vrtI][i];
 436     }
 437
 438     forAll(positions, posI)
 439     {
 440         neighbours[posI] = -1;
 441
 442         const label nei = findLeafLabelForPosition(positions[posI]);
 443
 444         if( (nei > -1) && leaves_[nei]->isLeaf() )
 445             neighbours[posI] = nei;
 446     }
 447
 448     # ifdef OCTREE_DEBUG
 449     Info << "Cube " << *oc << endl;
 450     Info << "Vertex " << short(vrtI) << endl;
 451     Info << "Neighbours " << endl;
 452     forAll(neighbours, nI)
 453         if( neighbours[nI] )
 454             Info << *neighbours[nI] << endl;
 455     # endif
 456 }
 457
 458 meshOctreeCube* meshOctree::findCubeForPosition
 459 (
 460     const meshOctreeCubeCoordinates& cc
 461 ) const
 462 {
 463     # ifdef OCTREE_DEBUG
 464     Info << "Finding position " << cc << endl;
 465     # endif
 466
 467     const label cpx = cc.posX();
 468     const label cpy = cc.posY();
 469     const label cpz = cc.posZ();
 470     const direction l = cc.level();
 471
 472     label levelLimiter = (1 << l);
 473     if(
 474         (cpx >= levelLimiter) || (cpx < 0) ||
 475         (cpy >= levelLimiter) || (cpy < 0) ||
 476         (!isQuadtree_ && (cpz >= levelLimiter || cpz < 0)) ||
 477         (isQuadtree_ && (cpz != initialCubePtr_->posZ()))
 478     )
 479     {
 480         return NULL;
 481     }
 482
 483     meshOctreeCube* neiPtr(initialCubePtr_);
 484
 485     for(label i=(l-1);i>=0;--i)
 486     {
 487         if( neiPtr && !neiPtr->isLeaf() )
 488         {
 489             levelLimiter = (1 << i);
 490
 491             label scI(0);
 492
 493             if( cpx & levelLimiter )
 494                 scI |= 1;
 495             if( cpy & levelLimiter )
 496                 scI |= 2;
 497             if( !isQuadtree_ && (cpz & levelLimiter) )
 498                 scI |= 4;
 499
 500             neiPtr = neiPtr->subCube(scI);
 501         }
 502         else
 503         {
 504             break;
 505         }
 506     }
 507
 508     # ifdef OCTREE_DEBUG
 509     Info << "Found position is " << *neiPtr << endl;
 510     # endif
 511
 512     return neiPtr;
 513 }
 514
 515 label meshOctree::findLeafLabelForPosition
 516 (
 517     const meshOctreeCubeCoordinates& cc
 518 ) const
 519 {
 520     const meshOctreeCube* ocPtr = findCubeForPosition(cc);
 521
 522     if( ocPtr && ocPtr->isLeaf() )
 523     {
 524         return ocPtr->cubeLabel();
 525     }
 526     else if( !ocPtr && (neiProcs_.size() != 0) )
 527     {
 528         const label levelLimiter = (1 << cc.level());
 529         if(
 530             (cc.posX() < levelLimiter) && (cc.posX() >= 0) &&
 531             (cc.posY() < levelLimiter) && (cc.posY() >= 0) &&
 532             ((!isQuadtree_ && (cc.posZ() < levelLimiter && cc.posZ() >= 0)) ||
 533             (isQuadtree_ && (cc.posZ() == initialCubePtr_->posZ())))
 534         )
 535         {
 536             return meshOctreeCubeBasic::OTHERPROC;
 537         }
 538     }
 539
 540     return -1;
 541 }
 542
 543 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 544
 545 } // End namespace Foam
 546
 547 // ************************************************************************* //