src/mesh/cfMesh/voronoiMesh/voronoiMeshExtractor/voronoiMeshExtractorCreateMesh.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::createPoints()
  45 {
  46     const LongList<point>& tetPoints = tetCreator_.tetPoints();
  47     const LongList<partTet>& tets = tetCreator_.tets();
  48
  49     pointFieldPMG& points = mesh_.points();
  50     points.setSize(tets.size());
  51
  52     # ifdef DEBUGVoronoi
  53     Info << "Number of tets " << tets.size() << endl;
  54     # endif
  55
  56     # ifdef USE_OMP
  57     # pragma omp parallel for schedule(dynamic, 100)
  58     # endif
  59     forAll(tets, tetI)
  60     {
  61         points[tetI] = tets[tetI].centroid(tetPoints);
  62
  63         # ifdef DEBUGVoronoi
  64         Info << "Point " << tetI << " has coordinates "
  65             << points[tetI] << endl;
  66         Info << "Tet of origin " << tetI << " has nodes "
  67             << tets[tetI] << endl;
  68         # endif
  69     }
  70 }
  71
  72 void voronoiMeshExtractor::createPolyMesh()
  73 {
  74     const VRWGraph& pointEdges = this->pointEdges();
  75     const VRWGraph& edgeTets = this->edgeTets();
  76     const boolList& boundaryEdge = this->boundaryEdge();
  77     const LongList<edge>& edges = this->edges();
  78
  79     polyMeshGenModifier meshModifier(mesh_);
  80     faceListPMG& faces = meshModifier.facesAccess();
  81     cellListPMG& cells = meshModifier.cellsAccess();
  82
  83     //- count the number of cells
  84     label nCells(0);
  85     labelList cellLabel(pointEdges.size(), -1);
  86
  87     forAll(pointEdges, pointI)
  88     {
  89         bool create(true);
  90         forAllRow(pointEdges, pointI, eI)
  91             if( boundaryEdge[pointEdges(pointI, eI)] )
  92             {
  93                 create = false;
  94                 break;
  95             }
  96
  97         if( !create || (pointEdges.sizeOfRow(pointI) == 0) )
  98             continue;
  99
 100         cellLabel[pointI] = nCells;
 101         ++nCells;
 102     }
 103
 104     # ifdef DEBUGVoronoi
 105     Info << "Number of cells " << nCells << endl;
 106     # endif
 107
 108     //- count the number of faces
 109     label nFaces(0);
 110     labelList faceLabel(edges.size(), -1);
 111
 112     forAll(boundaryEdge, edgeI)
 113     {
 114         if( boundaryEdge[edgeI] )
 115             continue;
 116
 117         const edge& e = edges[edgeI];
 118         if( cellLabel[e[0]] < 0 && cellLabel[e[1]] < 0 )
 119             continue;
 120
 121         faceLabel[edgeI] = nFaces;
 122         ++nFaces;
 123     }
 124
 125     # ifdef DEBUGVoronoi
 126     Info << "Number of mesh faces " << nFaces << endl;
 127     # endif
 128
 129     //- create faces
 130     Info << "Creating faces " << endl;
 131     faces.setSize(nFaces);
 132     labelList nFacesInCell(nCells, 0);
 133
 134     forAll(edges, edgeI)
 135     {
 136         if( faceLabel[edgeI] < 0 )
 137             continue;
 138
 139         const label faceI = faceLabel[edgeI];
 140
 141         face& f = faces[faceI];
 142         f.setSize(edgeTets.sizeOfRow(edgeI));
 143
 144         //- fill the faces with the node labels
 145         forAllRow(edgeTets, edgeI, pI)
 146             f[pI] = edgeTets(edgeI, pI);
 147
 148         const edge& e = edges[edgeI];
 149         const label cOwn = cellLabel[e[0]];
 150         const label cNei = cellLabel[e[1]];
 151
 152         if( cOwn < 0 || ((cOwn > cNei) && (cNei != -1)) )
 153         {
 154             # ifdef DEBUGVoronoi
 155             Info << "Reversing face " << cOwn << " " << cNei << endl;
 156             # endif
 157
 158             f = f.reverseFace();
 159         }
 160
 161         if( cOwn >= 0 )
 162             ++nFacesInCell[cOwn];
 163         if( cNei >= 0 )
 164             ++nFacesInCell[cNei];
 165     }
 166
 167     //- create cells
 168     # ifdef DEBUGVoronoi
 169     Info << "Setting cell sizes" << endl;
 170     # endif
 171
 172     cells.setSize(nCells);
 173     forAll(nFacesInCell, cellI)
 174     {
 175         cells[cellI].setSize(nFacesInCell[cellI]);
 176         nFacesInCell[cellI] = 0;
 177     }
 178
 179     # ifdef DEBUGVoronoi
 180     Info << "Filling cells" << endl;
 181     # endif
 182
 183     forAll(edges, edgeI)
 184     {
 185         if( faceLabel[edgeI] < 0 )
 186             continue;
 187
 188         const label faceI = faceLabel[edgeI];
 189         const edge& e = edges[edgeI];
 190         const label cOwn = cellLabel[e[0]];
 191         const label cNei = cellLabel[e[1]];
 192
 193         if( cOwn >= 0 )
 194             cells[cOwn][nFacesInCell[cOwn]++] = faceI;
 195         if( cNei >= 0 )
 196             cells[cNei][nFacesInCell[cNei]++] = faceI;
 197     }
 198
 199     # ifdef DEBUGVoronoi
 200     Info << "Finished generating cells" << endl;
 201     # endif
 202
 203     mesh_.clearAddressingData();
 204 }
 205
 206 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 207
 208 } // End namespace Foam
 209
 210 // ************************************************************************* //