src/mesh/cfMesh/tetMesh/tetMeshExtractorOctree/tetMeshExtractorOctree.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 "tetMeshExtractorOctree.H"
  29 #include "meshOctree.H"
  30 #include "triSurface.H"
  31 #include "polyMeshGenModifierAddCellByCell.H"
  32 #include "demandDrivenData.H"
  33
  34 # ifdef USE_OMP
  35 #include <omp.h>
  36 # endif
  37
  38 // #define DEBUGTets
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 void tetMeshExtractorOctree::createPoints()
  48 {
  49     polyMeshGenModifier meshModifier ( mesh_ );
  50     pointFieldPMG& points = meshModifier.pointsAccess();
  51
  52     const LongList<point>& tetPoints = tetCreator_.tetPoints();
  53
  54     points.setSize(tetPoints.size());
  55
  56     # ifdef USE_OMP
  57     # pragma omp parallel for
  58     # endif
  59     forAll(tetPoints, pointI)
  60         points[pointI] = tetPoints[pointI];
  61 }
  62
  63 void tetMeshExtractorOctree::createPolyMesh()
  64 {
  65     polyMeshGenModifier meshModifier ( mesh_ );
  66
  67     faceListPMG& faces = meshModifier.facesAccess();
  68     cellListPMG& cells = meshModifier.cellsAccess();
  69     meshModifier.boundariesAccess().setSize ( 0 );
  70     meshModifier.procBoundariesAccess().setSize ( 0 );
  71
  72     const LongList<partTet>& tets = tetCreator_.tets();
  73
  74     VRWGraph pTets;
  75     pTets.reverseAddressing(mesh_.points().size(), tets);
  76
  77     //- set the number of cells
  78     cells.setSize(tets.size());
  79
  80     //- all faces of tetrahedral cells
  81     faces.setSize(4*tets.size());
  82     boolList removeFace(faces.size());
  83
  84     # ifdef USE_OMP
  85     # pragma omp parallel if( tets.size() > 1000 )
  86     # endif
  87     {
  88         //- set face labels
  89         # ifdef USE_OMP
  90         # pragma omp for
  91         # endif
  92         forAll(removeFace, faceI)
  93             removeFace[faceI] = false;
  94
  95         //- set sizes of cells and create all faces
  96         # ifdef USE_OMP
  97         # pragma omp for schedule(dynamic, 20)
  98         # endif
  99         forAll(tets, elmtI)
 100         {
 101             cells[elmtI].setSize(4);
 102
 103             const partTet& elmt = tets[elmtI];
 104
 105             label faceI = 4 * elmtI;
 106
 107             //- first face
 108             cells[elmtI][0] = faceI;
 109
 110             face& f0 = faces[faceI];
 111             f0.setSize(3);
 112
 113             f0[0] = elmt.a();
 114             f0[1] = elmt.c();
 115             f0[2] = elmt.b();
 116
 117             ++faceI;
 118
 119             //- second face
 120             cells[elmtI][1] = faceI;
 121
 122             face& f1 = faces[faceI];
 123             f1.setSize(3);
 124
 125             f1[0] = elmt.a();
 126             f1[1] = elmt.b();
 127             f1[2] = elmt.d();
 128
 129             ++faceI;
 130
 131             //- third face
 132             cells[elmtI][2] = faceI;
 133
 134             face& f2 = faces[faceI];
 135             f2.setSize ( 3 );
 136
 137             f2[0] = elmt.b();
 138             f2[1] = elmt.c();
 139             f2[2] = elmt.d();
 140
 141             ++faceI;
 142
 143             //- fourth face
 144             cells[elmtI][3] = faceI;
 145
 146             face& f3 = faces[faceI];
 147             f3.setSize ( 3 );
 148
 149             f3[0] = elmt.c();
 150             f3[1] = elmt.a();
 151             f3[2] = elmt.d();
 152         }
 153
 154         # ifdef USE_OMP
 155         # pragma omp barrier
 156         # endif
 157
 158         //- find duplicate faces
 159         # ifdef USE_OMP
 160         # pragma omp for schedule(dynamic, 20)
 161         # endif
 162         forAll(cells, cellI)
 163         {
 164             cell& c = cells[cellI];
 165
 166             forAll(c, fI)
 167             {
 168                 const face& f = faces[c[fI]];
 169                 const label pointI = f[0];
 170
 171                 forAllRow(pTets, pointI, ptI)
 172                 {
 173                     //- do not check cells with greater labels
 174                     //- they cannot be face owners
 175                     if( pTets(pointI, ptI) >= cellI )
 176                         continue;
 177
 178                     const cell& otherTet = cells[pTets(pointI, ptI)];
 179
 180                     //- check faces created from a tet
 181                     forAll(otherTet, ofI)
 182                     {
 183                         //- do not compare faces with greater labels
 184                         //- they shall not be removed here
 185                         if( otherTet[ofI] >= c[fI] )
 186                             continue;
 187
 188                         //- check if the faces are equal
 189                         if( f == faces[otherTet[ofI]] )
 190                         {
 191                             removeFace[c[fI]] = true;
 192                             c[fI] = otherTet[ofI];
 193                         }
 194                     }
 195                 }
 196             }
 197         }
 198     }
 199
 200     //- remove duplicate faces
 201     label nFaces(0);
 202     labelLongList newFaceLabel(faces.size(), -1);
 203
 204     forAll(faces, faceI)
 205     {
 206         if( !removeFace[faceI] )
 207         {
 208             if( nFaces < faceI )
 209                 faces[nFaces].transfer(faces[faceI]);
 210
 211             newFaceLabel[faceI] = nFaces;
 212             ++nFaces;
 213         }
 214     }
 215
 216     //- set the size of faces
 217     faces.setSize(nFaces);
 218
 219     //- change cells
 220     # ifdef USE_OMP
 221     # pragma omp for schedule(dynamic, 40)
 222     # endif
 223     forAll(cells, cellI)
 224     {
 225         cell& c = cells[cellI];
 226
 227         DynList<label> newC;
 228
 229         forAll(c, fI)
 230         {
 231             if( newFaceLabel[c[fI]] != -1 )
 232                 newC.append(newFaceLabel[c[fI]]);
 233         }
 234
 235         c.setSize(newC.size());
 236         forAll(c, fI)
 237         c[fI] = newC[fI];
 238     }
 239 }
 240
 241 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 242
 243 // Construct from octree and mesh data
 244 tetMeshExtractorOctree::tetMeshExtractorOctree
 245 (
 246     const meshOctree& octree,
 247     const IOdictionary& meshDict,
 248     polyMeshGen& mesh
 249 )
 250 :
 251     tetCreator_(octree, meshDict),
 252     mesh_(mesh)
 253 {}
 254
 255 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 256
 257 tetMeshExtractorOctree::~tetMeshExtractorOctree()
 258 {}
 259
 260 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 261
 262 void tetMeshExtractorOctree::createMesh()
 263 {
 264     Info << "Extracting tetMesh" << endl;
 265
 266     //- copy tet points into the mesh
 267     createPoints();
 268
 269     //- create the mesh
 270     createPolyMesh();
 271
 272     polyMeshGenModifier(mesh_).reorderBoundaryFaces();
 273     polyMeshGenModifier(mesh_).removeUnusedVertices();
 274
 275     Info << "Mesh has :" << nl
 276     << mesh_.points().size() << " vertices " << nl
 277     << mesh_.faces().size() << " faces" << nl
 278     << mesh_.cells().size() << " cells" << endl;
 279
 280     Info << "Finished extracting tetMesh" << endl;
 281 }
 282
 283 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 284
 285 } // End namespace Foam
 286
 287 // ************************************************************************* //