src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.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 "demandDrivenData.H"
  29 #include "meshOptimizer.H"
  30 #include "meshSurfaceEngine.H"
  31 #include "meshSurfacePartitioner.H"
  32 #include "polyMeshGenAddressing.H"
  33 #include "polyMeshGenChecks.H"
  34 #include "labelLongList.H"
  35
  36 // #define DEBUGSmoothing
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * Private member functions * * * * * * * * * * * * * * * * * //
  44
  45 const meshSurfaceEngine& meshOptimizer::meshSurface() const
  46 {
  47     if( !msePtr_ )
  48         msePtr_ = new meshSurfaceEngine(mesh_);
  49
  50     return *msePtr_;
  51 }
  52
  53 void meshOptimizer::clearSurface()
  54 {
  55     deleteDemandDrivenData(msePtr_);
  56 }
  57
  58 label meshOptimizer::findBadFaces
  59 (
  60     labelHashSet& badFaces,
  61     const boolList& changedFace
  62 ) const
  63 {
  64     badFaces.clear();
  65
  66     polyMeshGenChecks::checkFacePyramids
  67     (
  68         mesh_,
  69         false,
  70         VSMALL,
  71         &badFaces,
  72         &changedFace
  73     );
  74
  75     polyMeshGenChecks::checkFaceFlatness
  76     (
  77         mesh_,
  78         false,
  79         0.8,
  80         &badFaces,
  81         &changedFace
  82     );
  83
  84     polyMeshGenChecks::checkCellPartTetrahedra
  85     (
  86         mesh_,
  87         false,
  88         VSMALL,
  89         &badFaces,
  90         &changedFace
  91     );
  92
  93     polyMeshGenChecks::checkFaceAreas
  94     (
  95         mesh_,
  96         false,
  97         VSMALL,
  98         &badFaces,
  99         &changedFace
 100     );
 101
 102     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
 103
 104     return nBadFaces;
 105 }
 106
 107 label meshOptimizer::findLowQualityFaces
 108 (
 109     labelHashSet& badFaces,
 110     const boolList& changedFace
 111 ) const
 112 {
 113     badFaces.clear();
 114
 115     polyMeshGenChecks::checkFaceDotProduct
 116     (
 117         mesh_,
 118         false,
 119         70.0,
 120         &badFaces
 121     );
 122
 123     polyMeshGenChecks::checkFaceSkewness
 124     (
 125         mesh_,
 126         false,
 127         2.0,
 128         &badFaces
 129     );
 130
 131     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
 132
 133     return nBadFaces;
 134 }
 135
 136 void meshOptimizer::calculatePointLocations()
 137 {
 138     vertexLocation_.setSize(mesh_.points().size());
 139     vertexLocation_ = INSIDE;
 140
 141     const meshSurfaceEngine& mse = meshSurface();
 142     const labelList& bPoints = mse.boundaryPoints();
 143
 144     //- mark boundary vertices
 145     forAll(bPoints, bpI)
 146         vertexLocation_[bPoints[bpI]] = BOUNDARY;
 147
 148     //- mark edge vertices
 149     meshSurfacePartitioner mPart(mse);
 150     forAllConstIter(labelHashSet, mPart.edgePoints(), it)
 151         vertexLocation_[bPoints[it.key()]] = EDGE;
 152
 153     //- mark corner vertices
 154     forAllConstIter(labelHashSet, mPart.corners(), it)
 155         vertexLocation_[bPoints[it.key()]] = CORNER;
 156
 157     if( Pstream::parRun() )
 158     {
 159         const polyMeshGenAddressing& addresing = mesh_.addressingData();
 160         const VRWGraph& pointAtProcs = addresing.pointAtProcs();
 161
 162         forAll(pointAtProcs, pointI)
 163             if( pointAtProcs.sizeOfRow(pointI) != 0 )
 164                 vertexLocation_[pointI] |= PARALLELBOUNDARY;
 165     }
 166 }
 167
 168 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 169
 170 // Construct from mesh
 171 meshOptimizer::meshOptimizer(polyMeshGen& mesh)
 172 :
 173     mesh_(mesh),
 174     vertexLocation_(),
 175     lockedFaces_(),
 176     msePtr_(NULL),
 177     enforceConstraints_(false),
 178     badPointsSubsetName_()
 179 {
 180     calculatePointLocations();
 181 }
 182
 183 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 184
 185 meshOptimizer::~meshOptimizer()
 186 {
 187     clearSurface();
 188 }
 189
 190 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 191
 192 void meshOptimizer::enforceConstraints(const word subsetName)
 193 {
 194     enforceConstraints_ = true;
 195
 196     badPointsSubsetName_ = subsetName;
 197 }
 198
 199 void meshOptimizer::lockCellsInSubset(const word& subsetName)
 200 {
 201     //- lock the points in the cell subset with the given name
 202     label subsetI = mesh_.cellSubsetIndex(subsetName);
 203     if( subsetI >= 0 )
 204     {
 205         labelLongList lc;
 206         mesh_.cellsInSubset(subsetI, lc);
 207
 208         lockCells(lc);
 209     }
 210     else
 211     {
 212         Warning << "Subset " << subsetName << " is not a cell subset!"
 213             << " Cannot lock cells!" << endl;
 214     }
 215 }
 216
 217 void meshOptimizer::lockFacesInSubset(const word& subsetName)
 218 {
 219     //- lock the points in the face subset with the given name
 220     label subsetI = mesh_.faceSubsetIndex(subsetName);
 221     if( subsetI >= 0 )
 222     {
 223         labelLongList lf;
 224         mesh_.facesInSubset(subsetI, lf);
 225
 226         lockFaces(lf);
 227     }
 228     else
 229     {
 230         Warning << "Subset " << subsetName << " is not a face subset!"
 231             << " Cannot lock faces!" << endl;
 232     }
 233 }
 234
 235 void meshOptimizer::lockPointsInSubset(const word& subsetName)
 236 {
 237     //- lock the points in the point subset with the given name
 238     label subsetI = mesh_.pointSubsetIndex(subsetName);
 239     if( subsetI >= 0 )
 240     {
 241         labelLongList lp;
 242         mesh_.pointsInSubset(subsetI, lp);
 243
 244         lockCells(lp);
 245     }
 246     else
 247     {
 248         Warning << "Subset " << subsetName << " is not a point subset!"
 249             << " Cannot lock points!" << endl;
 250     }
 251 }
 252
 253 void meshOptimizer::removeUserConstraints()
 254 {
 255     lockedFaces_.setSize(0);
 256
 257     //- unlock points
 258     # ifdef USE_OMP
 259     # pragma omp parallel for schedule(dynamic, 50)
 260     # endif
 261     forAll(vertexLocation_, i)
 262     {
 263         if( vertexLocation_[i] & LOCKED )
 264             vertexLocation_[i] ^= LOCKED;
 265     }
 266 }
 267
 268 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 269
 270 } // End namespace Foam
 271
 272 // ************************************************************************* //