src/mesh/cfMesh/utilities/intersectionTools/findCellsIntersectingSurface/findCellsIntersectingSurface.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 "findCellsIntersectingSurface.H"
  29 #include "polyMeshGen.H"
  30 #include "polyMeshGenAddressing.H"
  31 #include "triSurf.H"
  32 #include "boundBox.H"
  33 #include "helperFunctions.H"
  34 #include "meshOctree.H"
  35 #include "meshOctreeCreator.H"
  36 #include "HashSet.H"
  37
  38 # ifdef USE_OMP
  39 #include <omp.h>
  40 # endif
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 void findCellsIntersectingSurface::generateOctree(const triSurf& surf)
  50 {
  51     octreePtr_ = new meshOctree(surf);
  52
  53     meshOctreeCreator(*octreePtr_).createOctreeWithRefinedBoundary(15, 15);
  54 }
  55
  56 void findCellsIntersectingSurface::findIntersectedCells()
  57 {
  58     const pointFieldPMG& points = mesh_.points();
  59     const faceListPMG& faces = mesh_.faces();
  60     const cellListPMG& cells = mesh_.cells();
  61     const labelList& owner = mesh_.owner();
  62
  63     const vectorField& faceCentres = mesh_.addressingData().faceCentres();
  64     const vectorField& cellCentres = mesh_.addressingData().cellCentres();
  65
  66     meshOctreeModifier octreeModifier(*octreePtr_);
  67     const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
  68
  69     intersectedCells_.setSize(cells.size());
  70     facetsIntersectingCell_.setSize(cells.size());
  71
  72     const triSurf& surf = octreePtr_->surface();
  73     const VRWGraph& pointFacets = surf.pointFacets();
  74     const pointField& sp = surf.points();
  75
  76     # ifdef USE_OMP
  77     # pragma omp parallel for schedule(dynamic, 40)
  78     # endif
  79     forAll(cells, cellI)
  80     {
  81         bool intersected(false);
  82
  83         const cell& c = cells[cellI];
  84
  85         //- find the bounding box of the cell
  86         boundBox bb(cellCentres[cellI], cellCentres[cellI]);
  87
  88         forAll(c, fI)
  89         {
  90             const face& f = faces[c[fI]];
  91
  92             forAll(f, pI)
  93             {
  94                 bb.min() = Foam::min(bb.min(), points[f[pI]]);
  95                 bb.max() = Foam::max(bb.max(), points[f[pI]]);
  96             }
  97         }
  98
  99         const vector spanCorr = 0.01 * bb.span();
 100         bb.max() += spanCorr;
 101         bb.min() -= spanCorr;
 102
 103         //- find surface triangles within the bounding box
 104         DynList<label> leavesInBox;
 105         octreePtr_->findLeavesContainedInBox(bb, leavesInBox);
 106         labelHashSet triangles;
 107         forAll(leavesInBox, boxI)
 108         {
 109             const meshOctreeCube& oc = *leaves[leavesInBox[boxI]];
 110
 111             if( oc.hasContainedElements() )
 112             {
 113                 const meshOctreeSlot& slot = *oc.slotPtr();
 114                 const label ceI = oc.containedElements();
 115
 116                 forAllRow(slot.containedTriangles_, ceI, tI)
 117                     triangles.insert(slot.containedTriangles_(ceI, tI));
 118             }
 119         }
 120
 121         //- remove triangles which do not intersect the bounding box
 122         labelHashSet reasonableCandidates;
 123
 124         forAllConstIter(labelHashSet, triangles, tIter)
 125         {
 126             const labelledTri& tri = surf[tIter.key()];
 127
 128             boundBox obb(sp[tri[0]], sp[tri[0]]);
 129             for(label i=1;i<3;++i)
 130             {
 131                 const point& v = sp[tri[i]];
 132                 obb.min() = Foam::min(obb.min(), v);
 133                 obb.max() = Foam::max(obb.max(), v);
 134             }
 135
 136             const vector spanTriCorr = SMALL * obb.span();
 137             obb.min() -= spanTriCorr;
 138             obb.max() += spanTriCorr;
 139
 140             if( obb.overlaps(bb) )
 141                 reasonableCandidates.insert(tIter.key());
 142         }
 143
 144         triangles.transfer(reasonableCandidates);
 145
 146         //- check if any of the surface vertices is contained within the cell
 147         labelHashSet nodes, facetsInCell;
 148         forAllConstIter(labelHashSet, triangles, tIter)
 149         {
 150             const labelledTri& tri = surf[tIter.key()];
 151
 152             forAll(tri, i)
 153                 nodes.insert(tri[i]);
 154         }
 155
 156         //- check which surface nodes are within the cell
 157         forAllConstIter(labelHashSet, nodes, nIter)
 158         {
 159             const point& p = sp[nIter.key()];
 160
 161             if( !bb.contains(p) )
 162                 continue;
 163
 164             bool foundIntersection(false);
 165             forAll(c, fI)
 166             {
 167                 const face& f = faces[c[fI]];
 168
 169                 if( owner[c[fI]] == cellI )
 170                 {
 171                     forAll(f, pI)
 172                     {
 173                         tetrahedron<point, point> tet
 174                         (
 175                             points[f[pI]],
 176                             points[f.prevLabel(pI)],
 177                             faceCentres[c[fI]],
 178                             cellCentres[cellI]
 179                         );
 180
 181                         if( help::pointInTetrahedron(p, tet) )
 182                         {
 183                             intersected = true;
 184                             forAllRow(pointFacets, nIter.key(), ptI)
 185                             {
 186                                 facetsInCell.insert
 187                                 (
 188                                     pointFacets(nIter.key(), ptI)
 189                                 );
 190                             }
 191
 192                             foundIntersection = true;
 193                             break;
 194                         }
 195                     }
 196
 197                     if( foundIntersection )
 198                         break;
 199                 }
 200                 else
 201                 {
 202                     forAll(f, pI)
 203                     {
 204                         tetrahedron<point, point> tet
 205                         (
 206                             points[f[pI]],
 207                             points[f.nextLabel(pI)],
 208                             faceCentres[c[fI]],
 209                             cellCentres[cellI]
 210                         );
 211
 212                         if( help::pointInTetrahedron(p, tet) )
 213                         {
 214                             intersected = true;
 215                             forAllRow(pointFacets, nIter.key(), ptI)
 216                             {
 217                                 facetsInCell.insert
 218                                 (
 219                                     pointFacets(nIter.key(), ptI)
 220                                 );
 221                             }
 222
 223                             foundIntersection = true;
 224                             break;
 225                         }
 226                     }
 227
 228                     if( foundIntersection )
 229                         break;
 230                 }
 231             }
 232         }
 233
 234         //- check if any triangle in the surface mesh
 235         //- intersects any of the cell's faces
 236         forAllConstIter(labelHashSet, triangles, tIter)
 237         {
 238             if( facetsInCell.found(tIter.key()) )
 239                 continue;
 240
 241             forAll(c, fI)
 242             {
 243                 const face& f = faces[c[fI]];
 244
 245                 const bool intersect =
 246                     help::doFaceAndTriangleIntersect
 247                     (
 248                         surf,
 249                         tIter.key(),
 250                         f,
 251                         points
 252                     );
 253
 254                 if( intersect )
 255                 {
 256                     intersected = true;
 257                     facetsInCell.insert(tIter.key());
 258                     break;
 259                 }
 260             }
 261         }
 262
 263         //- store the results for this cell
 264         intersectedCells_[cellI] = intersected;
 265         if( intersected )
 266         {
 267             # ifdef USE_OMP
 268             # pragma omp critical
 269             # endif
 270             {
 271                 facetsIntersectingCell_.setRow(cellI, facetsInCell.toc());
 272             }
 273         }
 274     }
 275 }
 276
 277 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 278
 279 findCellsIntersectingSurface::findCellsIntersectingSurface
 280 (
 281     polyMeshGen& mesh,
 282     const meshOctree& octree
 283 )
 284 :
 285     mesh_(mesh),
 286     octreePtr_(const_cast<meshOctree*>(&octree)),
 287     octreeGenerated_(false),
 288     intersectedCells_(),
 289     facetsIntersectingCell_()
 290 {
 291     findIntersectedCells();
 292 }
 293
 294 findCellsIntersectingSurface::findCellsIntersectingSurface
 295 (
 296     polyMeshGen& mesh,
 297     const triSurf& surface
 298 )
 299 :
 300     mesh_(mesh),
 301     octreePtr_(NULL),
 302     octreeGenerated_(true),
 303     intersectedCells_(),
 304     facetsIntersectingCell_()
 305 {
 306     generateOctree(surface);
 307
 308     findIntersectedCells();
 309 }
 310
 311 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 312
 313 findCellsIntersectingSurface::~findCellsIntersectingSurface()
 314 {
 315     if( octreeGenerated_ )
 316         deleteDemandDrivenData(octreePtr_);
 317 }
 318
 319 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 320
 321 const boolList& findCellsIntersectingSurface::intersectedCells() const
 322 {
 323     return intersectedCells_;
 324 }
 325
 326 const VRWGraph& findCellsIntersectingSurface::facetsIntersectingCells() const
 327 {
 328     return facetsIntersectingCell_;
 329 }
 330
 331 void findCellsIntersectingSurface::addIntersectedCellsToSubset
 332 (
 333     const word subsetName
 334 )
 335 {
 336     const label setId = mesh_.addCellSubset(subsetName);
 337
 338     forAll(intersectedCells_, cellI)
 339         if( intersectedCells_[cellI] )
 340             mesh_.addCellToSubset(setId, cellI);
 341 }
 342
 343 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 344
 345 } // End namespace Foam
 346
 347 // ************************************************************************* //