src/meshLibrary/utilities/triSurfaceTools/triSurfacePartitioner/triSurfacePartitionerCreateAddressing.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 "triSurfacePartitioner.H"
  29 #include "demandDrivenData.H"
  30 #include "labelLongList.H"
  31 #include "boolList.H"
  32
  33 # ifdef USE_OMP
  34 #include <omp.h>
  35 # endif
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 void triSurfacePartitioner::calculatePatchAddressing()
  45 {
  46     calculateCornersAndAddressing();
  47
  48     calculatePatchPatches();
  49
  50     calculateEdgeGroups();
  51
  52     calculatePatchToEdgeGroups();
  53
  54     calculateEdgeGroupsToCorners();
  55 }
  56
  57 void triSurfacePartitioner::calculateCornersAndAddressing()
  58 {
  59     const VRWGraph& pointFaces = surface_.pointFacets();
  60     const edgeLongList& edges = surface_.edges();
  61     const VRWGraph& edgeFaces = surface_.edgeFacets();
  62
  63     //- find the number of feature edges connected to each surface node
  64     labelList nEdgesAtNode(surface_.points().size(), 0);
  65     forAll(edgeFaces, eI)
  66     {
  67         if( edgeFaces.sizeOfRow(eI) != 2 )
  68             continue;
  69
  70         const label sPatch = surface_[edgeFaces(eI, 0)].region();
  71         const label ePatch = surface_[edgeFaces(eI, 1)].region();
  72
  73         if( sPatch != ePatch )
  74         {
  75             const edge& e = edges[eI];
  76             ++nEdgesAtNode[e.start()];
  77             ++nEdgesAtNode[e.end()];
  78         }
  79     }
  80
  81     //- count the number of feature edges connected to each surface point
  82     //- corners must have 3 or more edges attached to them
  83     label nCorners(0);
  84     forAll(nEdgesAtNode, pI)
  85     {
  86         if( nEdgesAtNode[pI] < 3 )
  87             continue;
  88
  89         ++nCorners;
  90     }
  91
  92     corners_.setSize(nCorners);
  93     cornerPatches_.setSize(nCorners);
  94     nCorners = 0;
  95
  96     //- store corner data
  97     DynList<label> patches;
  98     forAll(pointFaces, pointI)
  99     {
 100         if( nEdgesAtNode[pointI] < 3 )
 101             continue;
 102
 103         patches.clear();
 104         forAllRow(pointFaces, pointI, pfI)
 105             patches.appendIfNotIn(surface_[pointFaces(pointI, pfI)].region());
 106
 107         corners_[nCorners] = pointI;
 108         cornerPatches_[nCorners] = patches;
 109         ++nCorners;
 110     }
 111 }
 112
 113 void triSurfacePartitioner::calculatePatchPatches()
 114 {
 115     const VRWGraph& edgeFaces = surface_.edgeFacets();
 116
 117     forAll(edgeFaces, eI)
 118     {
 119         if( edgeFaces.sizeOfRow(eI) != 2 )
 120         {
 121             Warning << "Surface is not a manifold!!" << endl;
 122             continue;
 123         }
 124
 125         const label sPatch = surface_[edgeFaces(eI, 0)].region();
 126         const label ePatch = surface_[edgeFaces(eI, 1)].region();
 127
 128         if( sPatch != ePatch )
 129         {
 130             patchPatches_[sPatch].insert(ePatch);
 131             patchPatches_[ePatch].insert(sPatch);
 132         }
 133     }
 134 }
 135
 136 void triSurfacePartitioner::calculateEdgeGroups()
 137 {
 138     const edgeLongList& edges = surface_.edges();
 139     const VRWGraph& edgeFaces = surface_.edgeFacets();
 140     const VRWGraph& pointEdges = surface_.pointEdges();
 141
 142
 143     //- make all feature edges
 144     boolList featureEdge(edgeFaces.size(), false);
 145
 146     # ifdef USE_OMP
 147     # pragma omp parallel for schedule(dynamic, 40)
 148     # endif
 149     forAll(edgeFaces, eI)
 150     {
 151         DynList<label> patches;
 152         forAllRow(edgeFaces, eI, efI)
 153             patches.appendIfNotIn(surface_[edgeFaces(eI, efI)].region());
 154
 155         if( patches.size() > 1 )
 156             featureEdge[eI] = true;
 157     }
 158
 159     //- create a set containing corners for fast searching
 160     labelHashSet corners;
 161     forAll(corners_, i)
 162         corners.insert(corners_[i]);
 163
 164     edgeGroups_.setSize(edgeFaces.size());
 165     edgeGroups_ = -1;
 166
 167     label nGroups(0);
 168     forAll(featureEdge, eI)
 169     {
 170         if( !featureEdge[eI] )
 171             continue;
 172         if( edgeGroups_[eI] != -1 )
 173             continue;
 174
 175         labelLongList front;
 176         front.append(eI);
 177         edgeGroups_[eI] = nGroups;
 178         featureEdge[eI] = false;
 179
 180         while( front.size() )
 181         {
 182             const label eLabel = front.removeLastElement();
 183             const edge& e = edges[eLabel];
 184
 185             for(label pI=0;pI<2;++pI)
 186             {
 187                 const label pointI = e[pI];
 188
 189                 if( corners.found(pointI) )
 190                     continue;
 191
 192                 forAllRow(pointEdges, pointI, peI)
 193                 {
 194                     const label eJ = pointEdges(pointI, peI);
 195
 196                     if( featureEdge[eJ] && (edgeGroups_[eJ] == -1) )
 197                     {
 198                         edgeGroups_[eJ] = nGroups;
 199                         featureEdge[eJ] = false;
 200                         front.append(eJ);
 201                     }
 202                 }
 203             }
 204         }
 205
 206         ++nGroups;
 207     }
 208
 209     Info << nGroups << " edge groups found!" << endl;
 210
 211     edgeGroupEdgeGroups_.clear();
 212     edgeGroupEdgeGroups_.setSize(nGroups);
 213 }
 214
 215 void triSurfacePartitioner::calculatePatchToEdgeGroups()
 216 {
 217     const VRWGraph& edgeFaces = surface_.edgeFacets();
 218
 219     forAll(edgeFaces, eI)
 220     {
 221         if( edgeGroups_[eI] < 0 )
 222             continue;
 223
 224         DynList<label> patches;
 225         forAllRow(edgeFaces, eI, efI)
 226             patches.appendIfNotIn(surface_[edgeFaces(eI, efI)].region());
 227
 228         forAll(patches, i)
 229         {
 230             const label patchI = patches[i];
 231
 232             for(label j=i+1;j<patches.size();++j)
 233             {
 234                 const label patchJ = patches[j];
 235
 236                 const std::pair<label, label> pp
 237                 (
 238                     Foam::min(patchI, patchJ),
 239                     Foam::max(patchI, patchJ)
 240                 );
 241
 242                 patchesEdgeGroups_[pp].insert(edgeGroups_[eI]);
 243             }
 244         }
 245     }
 246 }
 247
 248 void triSurfacePartitioner::calculateEdgeGroupsToCorners()
 249 {
 250     const VRWGraph& pointEdges = surface_.pointEdges();
 251
 252     forAll(corners_, cornerI)
 253     {
 254         DynList<label> edgeGroupsAtCorner;
 255         const label pointI = corners_[cornerI];
 256
 257         forAllRow(pointEdges, pointI, peI)
 258             edgeGroupsAtCorner.appendIfNotIn
 259             (
 260                 edgeGroups_[pointEdges(pointI, peI)]
 261             );
 262
 263         forAll(edgeGroupsAtCorner, i)
 264         {
 265             const label epI = edgeGroupsAtCorner[i];
 266             if( epI < 0 )
 267                 continue;
 268             for(label j=i+1;j<edgeGroupsAtCorner.size();++j)
 269             {
 270                 const label epJ = edgeGroupsAtCorner[j];
 271                 if( epJ < 0 )
 272                     continue;
 273
 274                 std::pair<label, label> ep
 275                 (
 276                     Foam::min(epI, epJ),
 277                     Foam::max(epI, epJ)
 278                 );
 279
 280                 //- create edgepartition - edge partitions addressing
 281                 edgeGroupEdgeGroups_[ep.first].insert(ep.second);
 282                 edgeGroupEdgeGroups_[ep.second].insert(ep.first);
 283
 284                 edgeGroupsCorners_[ep].insert(cornerI);
 285             }
 286         }
 287     }
 288 }
 289
 290 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 291
 292 } // End namespace Foam
 293
 294 // ************************************************************************* //