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