src/meshLibrary/utilities/triSurfaceTools/triSurfaceClassifyEdges/triSurfaceClassifyEdgesFunctions.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 "triSurfaceClassifyEdges.H"
  29 #include "demandDrivenData.H"
  30 #include "helperFunctions.H"
  31 #include "triSurf.H"
  32 #include "meshOctree.H"
  33 #include "labelPair.H"
  34
  35 #ifdef USE_OMP
  36 #include <omp.h>
  37 #endif
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 void triSurfaceClassifyEdges::checkOrientation()
  47 {
  48     const triSurf& surf = octree_.surface();
  49     const boundBox& rootBox = octree_.rootBox();
  50     const pointField& points = surf.points();
  51     const VRWGraph& facetEdges = surf.facetEdges();
  52     const VRWGraph& edgeFacets = surf.edgeFacets();
  53
  54     facetOrientation_.setSize(surf.size());
  55
  56     //- sort all surface facets into groups consisting of facets with consistent
  57     //- orientation. Do not cross non-manifold edges
  58     labelLongList orientationGroup(surf.size(), -1);
  59     label nGroups(0);
  60
  61     forAll(surf, triI)
  62     {
  63         if( orientationGroup[triI] != -1 )
  64             continue;
  65
  66         orientationGroup[triI] = nGroups;
  67         labelLongList front;
  68         front.append(triI);
  69
  70         while( front.size() != 0 )
  71         {
  72             const label tLabel = front.removeLastElement();
  73
  74             const labelledTri& facet = surf[tLabel];
  75
  76             forAll(facet, eI)
  77             {
  78                 const label edgeI = facetEdges(tLabel, eI);
  79
  80                 if( edgeFacets.sizeOfRow(edgeI) != 2 )
  81                     continue;
  82
  83                 forAllRow(edgeFacets, edgeI, efI)
  84                 {
  85                     const label neiFacetI = edgeFacets(edgeI, efI);
  86
  87                     if( orientationGroup[neiFacetI] != -1 )
  88                         continue;
  89                     if( neiFacetI == tLabel )
  90                         continue;
  91
  92                     const labelledTri& neiFacet = surf[neiFacetI];
  93
  94                     //- check the orientation of triangles at this edge
  95                     //- check the sign of the angle
  96                     //- if the orientation  is not consistent
  97                     DynList<labelPair, 2> sharedIndices;
  98                     forAll(facet, i)
  99                     {
 100                         forAll(neiFacet, j)
 101                         {
 102                             if( facet[i] == neiFacet[j] )
 103                                 sharedIndices.append(labelPair(i, j));
 104                         }
 105                     }
 106
 107                     if( sharedIndices.size() == 2 )
 108                     {
 109                         const labelPair& pair0 = sharedIndices[0];
 110                         const labelPair& pair1 = sharedIndices[1];
 111                         if( ((pair0.first() + 1) % 3) == pair1.first() )
 112                         {
 113                             if( (pair1.second() + 1) % 3 == pair0.second() )
 114                             {
 115                                 orientationGroup[neiFacetI] = nGroups;
 116                                 front.append(neiFacetI);
 117                             }
 118                         }
 119                         else
 120                         {
 121                             if( (pair0.second() + 1) % 3 == pair1.second() )
 122                             {
 123                                 orientationGroup[neiFacetI] = nGroups;
 124                                 front.append(neiFacetI);
 125                             }
 126                         }
 127                     }
 128                 }
 129             }
 130         }
 131
 132         ++nGroups;
 133     }
 134
 135     Info << "Found " << nGroups
 136          << " groups of triangles with consistent orientation" << endl;
 137
 138     //- find the octree leaves containing each triangle
 139     VRWGraph triangleInLeaves(surf.size());
 140     labelLongList ntl(surf.size(), 0);
 141
 142     DynList<label> helper;
 143     for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 144     {
 145         helper.clear();
 146         octree_.containedTriangles(leafI, helper);
 147
 148         forAll(helper, i)
 149             ++ntl[helper[i]];
 150     }
 151
 152     forAll(ntl, triI)
 153         triangleInLeaves.setRowSize(triI, ntl[triI]);
 154
 155     ntl = 0;
 156     for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 157     {
 158         helper.clear();
 159         octree_.containedTriangles(leafI, helper);
 160
 161         forAll(helper, i)
 162         {
 163             const label triI = helper[i];
 164
 165             triangleInLeaves(triI, ntl[triI]++) = leafI;
 166         }
 167     }
 168
 169     //- check the orientation of all facets in a group and collect their votes
 170     DynList<labelPair> groupVotes;
 171     groupVotes.setSize(nGroups);
 172     groupVotes = labelPair(0, 0);
 173
 174     # ifdef USE_OMP
 175     # pragma omp parallel if( surf.size() > 1000 ) private(helper)
 176     # endif
 177     {
 178         DynList<labelPair> localVotes;
 179         localVotes.setSize(nGroups);
 180         localVotes = labelPair(0, 0);
 181
 182         # ifdef USE_OMP
 183         # pragma omp for schedule(dynamic, 40)
 184         # endif
 185         forAll(orientationGroup, triI)
 186         {
 187             const labelledTri& tri = surf[triI];
 188             const point c = tri.centre(points);
 189             vector n = tri.normal(points);
 190             const scalar magN = mag(n);
 191
 192             if( magN < VSMALL )
 193                 continue;
 194
 195             n /= magN;
 196
 197             //- find the OUTSIDE octree cubes in the vicinity of the triangle
 198             //- and check the orientation of the triangle
 199             forAllRow(triangleInLeaves, triI, tlI)
 200             {
 201                 const label leafI = triangleInLeaves(triI, tlI);
 202
 203                 octree_.findAllLeafNeighbours(leafI, helper);
 204
 205                 forAll(helper, i)
 206                 {
 207                     const label leafJ = helper[i];
 208
 209                     const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafJ);
 210
 211                     if( oc.cubeType() & meshOctreeCubeBasic::OUTSIDE )
 212                     {
 213                         const scalar length = 3.0 * oc.size(rootBox);
 214
 215                         point pMin, pMax;
 216                         oc.cubeBox(rootBox, pMin, pMax);
 217
 218                         const boundBox bb(pMin, pMax);
 219
 220                         //- check whether the ray casted from
 221                         //- the centre of the triangle intersects the cube
 222                         const point endPos = c + length * n;
 223                         const point endNeg = c - length * n;
 224
 225                         if( help::boundBoxLineIntersection(c, endPos, bb) )
 226                         {
 227                             //- found an intersection in the positive direction
 228                             ++localVotes[orientationGroup[triI]].first();
 229                         }
 230                         else if( help::boundBoxLineIntersection(c, endNeg, bb) )
 231                         {
 232                             //- found an intersection in the negative direction
 233                             ++localVotes[orientationGroup[triI]].second();
 234                         }
 235                     }
 236                 }
 237             }
 238         }
 239
 240         # ifdef USE_OMP
 241         # pragma omp critical(grouping)
 242         # endif
 243         {
 244             forAll(localVotes, groupI)
 245             {
 246                 groupVotes[groupI].first() += localVotes[groupI].first();
 247                 groupVotes[groupI].second() += localVotes[groupI].second();
 248             }
 249         }
 250     }
 251
 252     Info << "Before determining of orientation" << endl;
 253
 254     //- determine whether a group is oriented outward or inward
 255     List<direction> outwardGroup(nGroups, direction(0));
 256
 257     forAll(groupVotes, groupI)
 258     {
 259         if( groupVotes[groupI].first() > groupVotes[groupI].second() )
 260         {
 261             outwardGroup[groupI] = 1;
 262         }
 263         else if( groupVotes[groupI].first() < groupVotes[groupI].second() )
 264         {
 265             outwardGroup[groupI] = 2;
 266         }
 267     }
 268
 269     Info << "Here setting orientation" << endl;
 270     //- Finally, set the orientation of the normal
 271     facetOrientation_.setSize(surf.size());
 272     forAll(facetOrientation_, triI)
 273         facetOrientation_[triI] = outwardGroup[orientationGroup[triI]];
 274 }
 275
 276 void triSurfaceClassifyEdges::classifyEdgesTypes()
 277 {
 278     const triSurf& surf = octree_.surface();
 279     const pointField& points = surf.points();
 280     const VRWGraph& edgeFacets = surf.edgeFacets();
 281     const edgeLongList& edges = surf.edges();
 282     const VRWGraph& pointEdges = surf.pointEdges();
 283     const edgeLongList& featureEdges = surf.featureEdges();
 284
 285     edgeTypes_.setSize(edgeFacets.size());
 286     edgeTypes_ = NONE;
 287
 288     //- set feature edge types
 289     # ifdef USE_OMP
 290     # pragma omp parallel for schedule(dynamic, 20)
 291     # endif
 292     forAll(featureEdges, feI)
 293     {
 294         const edge& e = featureEdges[feI];
 295
 296         forAllRow(pointEdges, e.start(), peI)
 297         {
 298             if( edges[pointEdges(e.start(), peI)] == e )
 299                 edgeTypes_[pointEdges(e.start(), peI)] |= FEATUREEDGE;
 300         }
 301     }
 302
 303     //- Finally, check the edges and assign flags
 304     # ifdef USE_OMP
 305     # pragma omp parallel for schedule(dynamic, 40)
 306     # endif
 307     forAll(edgeFacets, edgeI)
 308     {
 309         if( edgeFacets.sizeOfRow(edgeI) != 2 )
 310         {
 311             //- this is not a manifold edge
 312             edgeTypes_[edgeI] = NONE;
 313             continue;
 314         }
 315
 316         //- surface is a manifold at this edge
 317         const label f0 = edgeFacets(edgeI, 0);
 318         const label f1 = edgeFacets(edgeI, 1);
 319
 320         const labelledTri& tri0 = surf[f0];
 321         const labelledTri& tri1 = surf[f1];
 322
 323         if( tri0.region() != tri1.region() )
 324             edgeTypes_[edgeI] |= FEATUREEDGE;
 325
 326         label apexNode(-1);
 327         forAll(tri1, pI)
 328         {
 329             bool found(false);
 330             forAll(tri0, pJ)
 331             {
 332                 if( tri0[pJ] == tri1[pI] )
 333                 {
 334                     found = true;
 335                     break;
 336                 }
 337
 338                 if( found )
 339                     break;
 340             }
 341
 342             if( !found )
 343             {
 344                 apexNode = tri1[pI];
 345                 break;
 346             }
 347         }
 348
 349         const tetrahedron<point, point> tet
 350         (
 351             points[tri0[0]],
 352             points[tri0[1]],
 353             points[tri0[2]],
 354             points[apexNode]
 355         );
 356
 357         const scalar vol = tet.mag();
 358
 359         if( facetOrientation_[f0] == 1 )
 360         {
 361             //- facet is outward oriented
 362             if( vol < -VSMALL )
 363             {
 364                 //- this is a convex edge
 365                 edgeTypes_[edgeI] |= CONVEXEDGE;
 366             }
 367             else if( vol > VSMALL )
 368             {
 369                 //- this is a concave edge
 370                 edgeTypes_[edgeI] |= CONCAVEEDGE;
 371             }
 372             else
 373             {
 374                 edgeTypes_[edgeI] |= FLATSURFACEEDGE;
 375             }
 376         }
 377         else if( facetOrientation_[f0] == 2 )
 378         {
 379             //- facet is inward oriented
 380             if( vol > VSMALL )
 381             {
 382                 //- this is a convex edge
 383                 edgeTypes_[edgeI] |= CONVEXEDGE;
 384             }
 385             else if( vol < -VSMALL )
 386             {
 387                 //- this is a concave edge
 388                 edgeTypes_[edgeI] |= CONCAVEEDGE;
 389             }
 390             else
 391             {
 392                 edgeTypes_[edgeI] |= FLATSURFACEEDGE;
 393             }
 394         }
 395     }
 396 }
 397
 398 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 399
 400 } // End namespace Foam
 401
 402 // ************************************************************************* //