src/mesh/cfMesh/voronoiMesh/voronoiMeshExtractor/voronoiMeshExtractorAddressing.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 "voronoiMeshExtractor.H"
  29 #include "demandDrivenData.H"
  30
  31 # ifdef USE_OMP
  32 #include <omp.h>
  33 # endif
  34
  35 //#define DEBUGVoronoi
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 void voronoiMeshExtractor::createAddressing() const
  45 {
  46     if( pointEdgesPtr_ || edgeTetsPtr_ || boundaryEdgePtr_ || edgesPtr_ )
  47         return;
  48
  49     pointEdgesPtr_ = new VRWGraph(tetCreator_.tetPoints().size());
  50     VRWGraph& pointEdges = *pointEdgesPtr_;
  51
  52     edgeTetsPtr_ = new VRWGraph();
  53     VRWGraph& edgeTets = *edgeTetsPtr_;
  54
  55     boundaryEdgePtr_ = new boolList();
  56     boolList& boundaryEdges = *boundaryEdgePtr_;
  57
  58     edgesPtr_ = new LongList<edge>();
  59     LongList<edge>& edges = *edgesPtr_;
  60
  61     //- create edges and edgeTets
  62     const LongList<partTet>& tets = tetCreator_.tets();
  63
  64     VRWGraph pointTets;
  65     pointTets.reverseAddressing(tets);
  66
  67     forAll(tets, tetI)
  68     {
  69         const FixedList<edge, 6> tetEdges = tets[tetI].edges();
  70
  71         forAll(tetEdges, eI)
  72         {
  73             const edge& e = tetEdges[eI];
  74             const label start = e.start();
  75
  76             const row endTets = pointTets[e.end()];
  77
  78             bool store(true);
  79
  80             DynList<label> eTets;
  81             forAllRow(pointTets, start, ptI)
  82             {
  83                 const label tetJ = pointTets(start, ptI);
  84
  85                 if( !endTets.contains(tetJ) )
  86                     continue;
  87
  88                 if( tetJ < tetI )
  89                 {
  90                     store = false;
  91                     break;
  92                 }
  93
  94                 eTets.append(tetJ);
  95             }
  96
  97             if( store )
  98             {
  99                 edgeTets.appendList(eTets);
 100
 101                 edges.append(e);
 102             }
 103         }
 104     }
 105
 106     # ifdef DEBUGVoronoi
 107     Info << "Edge tets " << edgeTets << endl;
 108     # endif
 109
 110     //- calculate point-edges addressing
 111     pointEdges.reverseAddressing(edges);
 112
 113     //- sort edge-tets in circular order
 114     boundaryEdges.setSize(edgeTets.size());
 115     boundaryEdges = false;
 116
 117     # ifdef USE_OMP
 118     # pragma omp parallel for schedule(dynamic, 100)
 119     # endif
 120     forAll(edgeTets, edgeI)
 121     {
 122         const edge& e = edges[edgeI];
 123         row eTets = edgeTets[edgeI];
 124
 125         # ifdef DEBUGVoronoi
 126         Info << "Edge " << edgeI << " has points " << e << endl;
 127         # endif
 128
 129         forAll(eTets, tetI)
 130         {
 131             const partTet& pt = tets[eTets[tetI]];
 132
 133             # ifdef DEBUGVoronoi
 134             Info << "Checking tet " << eTets[tetI] << " points " << pt << endl;
 135             # endif
 136
 137             //- find the face shared with the neighbour
 138             const FixedList<edge, 6> tetEdges = pt.edges();
 139
 140             label searchPoint(-1);
 141             forAll(tetEdges, eI)
 142             {
 143                 if( tetEdges[eI] == e )
 144                 {
 145                     if( tetEdges[eI].start() == e.start() )
 146                     {
 147                         searchPoint = pt[sameOrientation_[eI]];
 148                     }
 149                     else
 150                     {
 151                         searchPoint = pt[oppositeOrientation_[eI]];
 152                     }
 153
 154                     break;
 155                 }
 156             }
 157
 158             if( searchPoint < 0 )
 159                 FatalErrorIn
 160                 (
 161                     "void voronoiMeshExtractor::createAddressing() const"
 162                 ) << " invalid search point " << abort(FatalError);
 163
 164             bool found(false);
 165             forAll(eTets, i)
 166             {
 167                 if( tetI == i )
 168                     continue;
 169
 170                 const partTet& ptNei = tets[eTets[i]];
 171
 172                 const label pos = ptNei.whichPosition(searchPoint);
 173
 174                 if( pos < 0 )
 175                     continue;
 176
 177                 if( tetI < eTets.size()-1 )
 178                 {
 179                     const label add  = eTets[tetI+1];
 180                     eTets[tetI+1] = eTets[i];
 181                     eTets[i] = add;
 182                 }
 183                 found = true;
 184                 break;
 185             }
 186
 187             if( !found )
 188                 boundaryEdges[edgeI] = true;
 189         }
 190     }
 191 }
 192
 193 const VRWGraph& voronoiMeshExtractor::pointEdges() const
 194 {
 195     if( !pointEdgesPtr_ )
 196         createAddressing();
 197
 198     return *pointEdgesPtr_;
 199 }
 200
 201 const LongList<edge>& voronoiMeshExtractor::edges() const
 202 {
 203     if( !edgesPtr_ )
 204         createAddressing();
 205
 206     return *edgesPtr_;
 207 }
 208
 209 const VRWGraph& voronoiMeshExtractor::edgeTets() const
 210 {
 211     if( !edgeTetsPtr_ )
 212         createAddressing();
 213
 214     return *edgeTetsPtr_;
 215 }
 216
 217 const boolList& voronoiMeshExtractor::boundaryEdge() const
 218 {
 219     if( !boundaryEdgePtr_ )
 220         createAddressing();
 221
 222     return *boundaryEdgePtr_;
 223 }
 224
 225 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 226
 227 } // End namespace Foam
 228
 229 // ************************************************************************* //