src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAddressing/meshOctreeAddressing.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 "meshOctreeAddressing.H"
  29 #include "demandDrivenData.H"
  30 #include "IOdictionary.H"
  31 #include "helperFunctions.H"
  32 #include "triSurf.H"
  33 #include "meshOctreeModifier.H"
  34
  35 // #define DEBUGSearch
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 void meshOctreeAddressing::clearOut()
  43 {
  44     clearNodeAddressing();
  45     clearBoxTypes();
  46     clearOctreeFaces();
  47     clearAddressing();
  48     clearParallelAddressing();
  49 }
  50
  51 void meshOctreeAddressing::clearNodeAddressing()
  52 {
  53     nNodes_ = 0;
  54     deleteDemandDrivenData(octreePointsPtr_);
  55     deleteDemandDrivenData(nodeLabelsPtr_);
  56     deleteDemandDrivenData(nodeLeavesPtr_);
  57
  58     deleteDemandDrivenData(nodeTypePtr_);
  59 }
  60
  61 void meshOctreeAddressing::clearBoxTypes()
  62 {
  63     deleteDemandDrivenData(boxTypePtr_);
  64 }
  65
  66 void meshOctreeAddressing::clearOctreeFaces()
  67 {
  68     deleteDemandDrivenData(octreeFacesPtr_);
  69     deleteDemandDrivenData(octreeFacesOwnersPtr_);
  70     deleteDemandDrivenData(octreeFacesNeighboursPtr_);
  71 }
  72
  73 void meshOctreeAddressing::clearAddressing()
  74 {
  75     deleteDemandDrivenData(leafFacesPtr_);
  76     deleteDemandDrivenData(nodeFacesPtr_);
  77     deleteDemandDrivenData(leafLeavesPtr_);
  78     deleteDemandDrivenData(octreeEdgesPtr_);
  79     deleteDemandDrivenData(edgeLeavesPtr_);
  80     deleteDemandDrivenData(leafEdgesPtr_);
  81     deleteDemandDrivenData(nodeEdgesPtr_);
  82     deleteDemandDrivenData(faceEdgesPtr_);
  83     deleteDemandDrivenData(edgeFacesPtr_);
  84 }
  85
  86 void meshOctreeAddressing::clearParallelAddressing()
  87 {
  88     deleteDemandDrivenData(globalPointLabelPtr_);
  89     deleteDemandDrivenData(globalPointToLocalPtr_);
  90     deleteDemandDrivenData(pointProcsPtr_);
  91     deleteDemandDrivenData(globalFaceLabelPtr_);
  92     deleteDemandDrivenData(globalFaceToLocalPtr_);
  93     deleteDemandDrivenData(faceProcsPtr_);
  94     deleteDemandDrivenData(globalLeafLabelPtr_);
  95     deleteDemandDrivenData(globalLeafToLocalPtr_);
  96     deleteDemandDrivenData(leafAtProcsPtr_);
  97 }
  98
  99 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 100
 101 // Construct from octree and IOdictionary
 102 meshOctreeAddressing::meshOctreeAddressing
 103 (
 104     const meshOctree& mo,
 105     const dictionary& dict,
 106     bool useDATABoxes
 107 )
 108 :
 109     octree_(mo),
 110     meshDict_(dict),
 111     useDATABoxes_(useDATABoxes),
 112     nNodes_(0),
 113     octreePointsPtr_(NULL),
 114     nodeLabelsPtr_(NULL),
 115     nodeLeavesPtr_(NULL),
 116     boxTypePtr_(NULL),
 117     nodeTypePtr_(NULL),
 118     octreeFacesPtr_(NULL),
 119     octreeFacesOwnersPtr_(NULL),
 120     octreeFacesNeighboursPtr_(NULL),
 121     leafFacesPtr_(NULL),
 122     nodeFacesPtr_(NULL),
 123     leafLeavesPtr_(NULL),
 124     octreeEdgesPtr_(NULL),
 125     edgeLeavesPtr_(NULL),
 126     leafEdgesPtr_(NULL),
 127     nodeEdgesPtr_(NULL),
 128     faceEdgesPtr_(NULL),
 129     edgeFacesPtr_(NULL),
 130     globalPointLabelPtr_(NULL),
 131     globalPointToLocalPtr_(NULL),
 132     pointProcsPtr_(NULL),
 133     globalFaceLabelPtr_(NULL),
 134     globalFaceToLocalPtr_(NULL),
 135     faceProcsPtr_(NULL),
 136     globalLeafLabelPtr_(NULL),
 137     globalLeafToLocalPtr_(NULL),
 138     leafAtProcsPtr_(NULL)
 139 {
 140     if( !useDATABoxes && dict.found("keepCellsIntersectingBoundary") )
 141     {
 142         useDATABoxes_ = readBool(dict.lookup("keepCellsIntersectingBoundary"));
 143     }
 144
 145     if( dict.found("nonManifoldMeshing") )
 146     {
 147         const bool nonManifoldMesh
 148         (
 149             readBool(dict.lookup("nonManifoldMeshing"))
 150         );
 151
 152         if( nonManifoldMesh )
 153             useDATABoxes_ = true;
 154     }
 155
 156     if( Pstream::parRun() )
 157     {
 158         meshOctreeModifier om(const_cast<meshOctree&>(octree_));
 159         om.addLayerFromNeighbouringProcessors();
 160     }
 161
 162     //- check for glued regions
 163     checkGluedRegions();
 164 }
 165
 166 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 167
 168 meshOctreeAddressing::~meshOctreeAddressing()
 169 {
 170     clearOut();
 171 }
 172
 173 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 174
 175 bool meshOctreeAddressing::isIntersectedFace(const label fI) const
 176 {
 177     const labelLongList& owner = octreeFaceOwner();
 178     const labelLongList& neighbour = octreeFaceNeighbour();
 179
 180     if( neighbour[fI] < 0 )
 181         return false;
 182
 183     Map<label> nAppearances;
 184     DynList<label> triangles;
 185     octree_.containedTriangles(owner[fI], triangles);
 186     forAll(triangles, triI)
 187     {
 188         if( nAppearances.found(triangles[triI]) )
 189         {
 190             ++nAppearances[triangles[triI]];
 191         }
 192         else
 193         {
 194             nAppearances.insert(triangles[triI], 1);
 195         }
 196     }
 197
 198     triangles.clear();
 199     octree_.containedTriangles(neighbour[fI], triangles);
 200     forAll(triangles, triI)
 201     {
 202         if( nAppearances.found(triangles[triI]) )
 203         {
 204             ++nAppearances[triangles[triI]];
 205         }
 206         else
 207         {
 208             nAppearances.insert(triangles[triI], 1);
 209         }
 210     }
 211
 212     forAllConstIter(Map<label>, nAppearances, iter)
 213     {
 214         if( iter() == 2 )
 215         {
 216             if
 217             (
 218                 octree_.returnLeaf(owner[fI]).level() ==
 219                 octree_.returnLeaf(neighbour[fI]).level()
 220             )
 221                 return true;
 222
 223             //- check intersection by geometric testing
 224             const triSurf& surf = octree_.surface();
 225             const pointField& points = this->octreePoints();
 226             const VRWGraph& faces = this->octreeFaces();
 227
 228             face f(faces.sizeOfRow(fI));
 229             forAll(f, pI)
 230                 f[pI] = faces(fI, pI);
 231
 232             if(
 233                 help::doFaceAndTriangleIntersect
 234                 (
 235                     surf,
 236                     iter.key(),
 237                     f,
 238                     points
 239                 )
 240             )
 241                 return true;
 242         }
 243     }
 244
 245     return false;
 246 }
 247
 248 bool meshOctreeAddressing::isIntersectedEdge(const label eI) const
 249 {
 250     const VRWGraph& edgeCubes = this->edgeLeaves();
 251
 252     Map<label> nAppearances;
 253     DynList<label> triangles;
 254     bool sameLevel(true);
 255
 256     forAllRow(edgeCubes, eI, i)
 257     {
 258         const label leafI = edgeCubes(eI, i);
 259         if( !octree_.hasContainedTriangles(leafI) )
 260             return false;
 261
 262         if
 263         (
 264             octree_.returnLeaf(leafI).level() !=
 265             octree_.returnLeaf(edgeCubes(eI, 0)).level()
 266         )
 267             sameLevel = false;
 268
 269         triangles.clear();
 270         octree_.containedTriangles(leafI, triangles);
 271         forAll(triangles, triI)
 272         {
 273             if( nAppearances.found(triangles[triI]) )
 274             {
 275                 ++nAppearances[triangles[triI]];
 276             }
 277             else
 278             {
 279                 nAppearances.insert(triangles[triI], 1);
 280             }
 281         }
 282     }
 283
 284     forAllConstIter(Map<label>, nAppearances, iter)
 285     {
 286         if( iter() == edgeCubes.sizeOfRow(eI) )
 287         {
 288             if( sameLevel )
 289                 return true;
 290
 291             //- check for geometric intersection
 292             const LongList<edge>& edges = this->octreeEdges();
 293             const pointField& points = this->octreePoints();
 294             point intersection(vector::zero);
 295
 296             if(
 297                 help::triLineIntersection
 298                 (
 299                     octree_.surface(),
 300                     iter.key(),
 301                     points[edges[eI].start()],
 302                     points[edges[eI].end()],
 303                     intersection
 304                 )
 305             )
 306                 return true;
 307         }
 308     }
 309
 310     return false;
 311 }
 312
 313 void meshOctreeAddressing::edgeIntersections
 314 (
 315     const label eI,
 316     DynList<point>& intersections
 317 ) const
 318 {
 319     intersections.clear();
 320
 321     const LongList<edge>& edges = this->octreeEdges();
 322     const pointField& points = this->octreePoints();
 323     const VRWGraph& edgeCubes = this->edgeLeaves();
 324     const scalar tol =
 325         SMALL * mag(points[edges[eI].start()] - points[edges[eI].end()]);
 326
 327     Map<label> nAppearances;
 328     DynList<label> triangles;
 329
 330     forAllRow(edgeCubes, eI, i)
 331     {
 332         const label leafI = edgeCubes(eI, i);
 333         if( !octree_.hasContainedTriangles(leafI) )
 334             return;
 335
 336         triangles.clear();
 337         octree_.containedTriangles(leafI, triangles);
 338         forAll(triangles, triI)
 339         {
 340             if( nAppearances.found(triangles[triI]) )
 341             {
 342                 ++nAppearances[triangles[triI]];
 343             }
 344             else
 345             {
 346                 nAppearances.insert(triangles[triI], 1);
 347             }
 348         }
 349     }
 350
 351     point intersection(vector::zero);
 352
 353     forAllConstIter(Map<label>, nAppearances, iter)
 354     {
 355         if( iter() == edgeCubes.sizeOfRow(eI) )
 356         {
 357             //- check for geometric intersection
 358             const bool intersectionExists =
 359                 help::triLineIntersection
 360                 (
 361                     octree_.surface(),
 362                     iter.key(),
 363                     points[edges[eI].start()],
 364                     points[edges[eI].end()],
 365                     intersection
 366                 );
 367
 368             if( intersectionExists )
 369             {
 370                 bool store(true);
 371                 forAll(intersections, i)
 372                     if( mag(intersections[i] - intersection) <= tol )
 373                         store = false;
 374
 375                 if( store )
 376                     intersections.append(intersection);
 377             }
 378         }
 379     }
 380 }
 381
 382 void meshOctreeAddressing::cubesAroundEdge
 383 (
 384     const label leafI,
 385     const direction eI,
 386     FixedList<label, 4>& edgeCubes
 387 ) const
 388 {
 389     const VRWGraph& nl = this->nodeLabels();
 390     const label nodeI = nl(leafI, meshOctreeCubeCoordinates::edgeNodes_[eI][0]);
 391     const FRWGraph<label, 8>& pLeaves = this->nodeLeaves();
 392
 393     switch( eI )
 394     {
 395         case 0: case 1: case 2: case 3:
 396         {
 397             edgeCubes[0] = pLeaves(nodeI, 1);
 398             edgeCubes[1] = pLeaves(nodeI, 3);
 399             edgeCubes[2] = pLeaves(nodeI, 5);
 400             edgeCubes[3] = pLeaves(nodeI, 7);
 401         } break;
 402         case 4: case 5: case 6: case 7:
 403         {
 404             edgeCubes[0] = pLeaves(nodeI, 2);
 405             edgeCubes[1] = pLeaves(nodeI, 3);
 406             edgeCubes[2] = pLeaves(nodeI, 6);
 407             edgeCubes[3] = pLeaves(nodeI, 7);
 408         } break;
 409         case 8: case 9: case 10: case 11:
 410         {
 411             edgeCubes[0] = pLeaves(nodeI, 4);
 412             edgeCubes[1] = pLeaves(nodeI, 5);
 413             edgeCubes[2] = pLeaves(nodeI, 6);
 414             edgeCubes[3] = pLeaves(nodeI, 7);
 415         } break;
 416         default:
 417         {
 418             FatalErrorIn
 419             (
 420                 "void tetMeshExtractorOctree::cubesAroundEdge(const label,"
 421                 "const direction, FixedList<label, 4>&)"
 422             ) << "Invalid edge specified!!" << abort(FatalError);
 423         } break;
 424     };
 425 }
 426
 427 label meshOctreeAddressing::findEdgeCentre
 428 (
 429     const label leafI,
 430     const direction eI
 431 ) const
 432 {
 433     if( octree_.isQuadtree() && eI >= 8 )
 434         return -1;
 435
 436     const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
 437     const VRWGraph& nl = this->nodeLabels();
 438     const label nodeI = nl(leafI, meshOctreeCubeCoordinates::edgeNodes_[eI][0]);
 439     const FRWGraph<label, 8>& pLeaves = this->nodeLeaves();
 440
 441     const direction level = oc.level();
 442
 443     label fI(-1);
 444     switch( eI )
 445     {
 446         case 0: case 1: case 2: case 3:
 447         {
 448             fI = 1;
 449         } break;
 450         case 4: case 5: case 6: case 7:
 451         {
 452             fI = 3;
 453         } break;
 454         case 8: case 9: case 10: case 11:
 455         {
 456             fI = 5;
 457         } break;
 458         default:
 459         {
 460             FatalErrorIn
 461             (
 462                 "label meshOctreeAddressing::findEdgeCentre"
 463                 "(const label leafI, const direction eI) const"
 464             ) << "Invalid edge specified!!" << abort(FatalError);
 465         } break;
 466     };
 467
 468     for(label i=0;i<4;++i)
 469     {
 470         const label fNode = meshOctreeCubeCoordinates::faceNodes_[fI][i];
 471
 472         if( pLeaves(nodeI, fNode) < 0 )
 473             continue;
 474
 475         const label leafJ = pLeaves(nodeI, fNode);
 476         if( octree_.returnLeaf(leafJ).level() > level )
 477         {
 478             const label shift = (i+2)%4;
 479             return nl(leafJ, meshOctreeCubeCoordinates::faceNodes_[fI][shift]);
 480         }
 481     }
 482
 483     return -1;
 484 }
 485
 486 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 487
 488 } // End namespace Foam
 489
 490 // ************************************************************************* //