src/mesh/cfMesh/meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicates/triSurfaceCleanupDuplicatesFunctions.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 "triSurfaceCleanupDuplicates.H"
  29 #include "triSurfModifier.H"
  30 #include "meshOctree.H"
  31 #include "demandDrivenData.H"
  32
  33 # ifdef USE_OMP
  34 #include <omp.h>
  35 # endif
  36 #include <set>
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 bool triSurfaceCleanupDuplicates::checkDuplicateTriangles()
  46 {
  47     labelLongList newTriangleLabel(surf_.size(), -1);
  48
  49     const VRWGraph& pointTriangles = surf_.pointFacets();
  50
  51     //- check if there exist duplicate triangles
  52     label counter(0);
  53
  54     forAll(surf_, triI)
  55     {
  56         if( newTriangleLabel[triI] != -1 )
  57             continue;
  58
  59         newTriangleLabel[triI] = counter;
  60         ++counter;
  61
  62         const labelledTri& tri = surf_[triI];
  63
  64         forAll(pointTriangles[tri[0]], ptI)
  65         {
  66             const label triJ = pointTriangles(tri[0], ptI);
  67
  68             if( triJ <= triI )
  69                 continue;
  70
  71             const labelledTri& otherTri = surf_[triJ];
  72
  73             if( tri == otherTri )
  74                 newTriangleLabel[triJ] = newTriangleLabel[triI];
  75         }
  76     }
  77
  78     Info << "Found " << (newTriangleLabel.size()-counter)
  79         << " duplicate triangles" << endl;
  80
  81     //- return if there exist no duplicate triangles
  82     if( counter == newTriangleLabel.size() )
  83         return false;
  84
  85     Info << "Current number of triangles" << surf_.size() << endl;
  86     Info << "New number of triangles " << counter << endl;
  87
  88     //- create new list of triangles and store it in the surface mesh
  89     LongList<labelledTri> newTriangles(counter);
  90
  91     forAll(newTriangleLabel, triI)
  92     {
  93         newTriangles[newTriangleLabel[triI]] = surf_[triI];
  94     }
  95
  96     updateTriangleLabels(newTriangleLabel);
  97
  98     triSurfModifier(surf_).facetsAccess().transfer(newTriangles);
  99     surf_.updateFacetsSubsets(newTriangleLabel);
 100
 101     return true;
 102 }
 103
 104 bool triSurfaceCleanupDuplicates::mergeDuplicatePoints()
 105 {
 106     pointField& pts = const_cast<pointField&>(surf_.points());
 107     labelLongList newPointLabel(surf_.nPoints());
 108     bool foundDuplicates(false);
 109
 110     # ifdef USE_OMP
 111     # pragma omp parallel
 112     # endif
 113     {
 114         # ifdef USE_OMP
 115         # pragma omp for
 116         # endif
 117         forAll(newPointLabel, pI)
 118             newPointLabel[pI] = pI;
 119
 120         # ifdef USE_OMP
 121         # pragma omp barrier
 122         # endif
 123         //- check if there exist any vertices closer
 124         //- than the prescribed tolerance
 125         # ifdef USE_OMP
 126         # pragma omp for schedule(dynamic, 20)
 127         # endif
 128         for(label leafI=0;leafI<octree_.numberOfLeaves();++leafI)
 129         {
 130             DynList<label> ct;
 131             octree_.containedTriangles(leafI, ct);
 132
 133             std::set<label> points;
 134
 135             forAll(ct, ctI)
 136             {
 137                 const label triI = newTriangleLabel_[ct[ctI]];
 138
 139                 if( triI < 0 )
 140                     continue;
 141
 142                 const labelledTri& tri = surf_[triI];
 143
 144                 forAll(tri, i)
 145                     points.insert(tri[i]);
 146             }
 147
 148             for
 149             (
 150                 std::set<label>::const_iterator it=points.begin();
 151                 it!=points.end();
 152             )
 153             {
 154                 const label pointI = *it;
 155
 156                 for
 157                 (
 158                     std::set<label>::const_iterator nIt=++it;
 159                     nIt!=points.end();
 160                     ++nIt
 161                 )
 162                     if( magSqr(pts[pointI] - pts[*nIt]) < sqr(tolerance_) )
 163                     {
 164                         foundDuplicates = true;
 165                         # ifdef USE_OMP
 166                         # pragma omp critical
 167                         # endif
 168                         newPointLabel[*nIt] = pointI;
 169                     }
 170             }
 171         }
 172     }
 173
 174     //- find if there exist no duplicate points
 175     if( !foundDuplicates )
 176         return false;
 177
 178     //- remove vertices and update node labels
 179     label counter(0);
 180     forAll(pts, pI)
 181         if( newPointLabel[pI] == pI )
 182         {
 183             newPointLabel[pI] = counter;
 184             if( counter < pI )
 185                 pts[counter] = pts[pI];
 186             ++counter;
 187         }
 188         else
 189         {
 190             const label origI = newPointLabel[pI];
 191             newPointLabel[pI] = newPointLabel[origI];
 192         }
 193
 194     Info << "Found " << (pts.size() - counter) << "duplicate points" << endl;
 195
 196     pts.setSize(counter);
 197
 198     //- remove triangles containing duplicate points
 199     LongList<labelledTri> newTriangles(surf_.facets());
 200     labelLongList newTriangleLabel(surf_.size(), -1);
 201
 202     counter = 0;
 203     forAll(surf_, triI)
 204     {
 205         const labelledTri& tri = surf_[triI];
 206         const labelledTri newTri
 207         (
 208             newPointLabel[tri[0]],
 209             newPointLabel[tri[1]],
 210             newPointLabel[tri[2]],
 211             tri.region()
 212         );
 213
 214         bool store(true);
 215         for(label i=0;i<2;++i)
 216             for(label j=i+1;j<3;++j)
 217                 if( newTri[i] == newTri[j] )
 218                 {
 219                     store = false;
 220                     break;
 221                 }
 222
 223         if( store )
 224         {
 225             newTriangles[counter] = newTri;
 226             newTriangleLabel[triI] = counter;
 227             ++counter;
 228         }
 229     }
 230
 231     newTriangles.setSize(counter);
 232
 233     //- update the surface
 234     triSurfModifier(surf_).facetsAccess().transfer(newTriangles);
 235     surf_.updateFacetsSubsets(newTriangleLabel);
 236
 237     return true;
 238 }
 239
 240 void triSurfaceCleanupDuplicates::updateTriangleLabels
 241 (
 242     const labelLongList& newTriangleLabel
 243 )
 244 {
 245     //- update addressing between the original triangles and the cleaned mesh
 246     forAll(newTriangleLabel_, triI)
 247     {
 248         const label origI = newTriangleLabel_[triI];
 249
 250         if( origI >= 0 )
 251             newTriangleLabel_[triI] = newTriangleLabel[origI];
 252     }
 253 }
 254
 255 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 256
 257 } // End namespace Foam
 258
 259 // ************************************************************************* //