src/meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.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 "polyMeshGenAddressing.H"
  31 #include "partTetMesh.H"
  32 #include "HashSet.H"
  33
  34 #include "tetMeshOptimisation.H"
  35
  36 //#define DEBUGSmooth
  37
  38 # ifdef DEBUGSmooth
  39 #include "helperFunctions.H"
  40 #include "polyMeshGenModifier.H"
  41 # endif
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 void meshOptimizer::untangleMeshFV()
  51 {
  52     Info << "Starting untangling the mesh" << endl;
  53
  54     # ifdef DEBUGSmooth
  55     partTetMesh tm(mesh_);
  56     forAll(tm.tets(), tetI)
  57         if( tm.tets()[tetI].mag(tm.points()) < 0.0 )
  58             Info << "Tet " << tetI << " is inverted!" << endl;
  59     polyMeshGen tetPolyMesh(mesh_.returnTime());
  60     tm.createPolyMesh(tetPolyMesh);
  61     polyMeshGenModifier(tetPolyMesh).removeUnusedVertices();
  62     forAll(tm.smoothVertex(), pI)
  63         if( !tm.smoothVertex()[pI] )
  64             Info << "Point " << pI << " cannot be moved!" << endl;
  65
  66     const VRWGraph& pTets = tm.pointTets();
  67     forAll(pTets, pointI)
  68     {
  69         const LongList<partTet>& tets = tm.tets();
  70         forAllRow(pTets, pointI, i)
  71             if( tets[pTets(pointI, i)].whichPosition(pointI) < 0 )
  72                 FatalError << "Wrong partTet" << abort(FatalError);
  73
  74         partTetMeshSimplex simplex(tm, pointI);
  75     }
  76
  77     boolList boundaryVertex(tetPolyMesh.points().size(), false);
  78     const labelList& neighbour = tetPolyMesh.neighbour();
  79     forAll(neighbour, faceI)
  80         if( neighbour[faceI] == -1 )
  81         {
  82             const face& f = tetPolyMesh.faces()[faceI];
  83
  84             forAll(f, pI)
  85                 boundaryVertex[f[pI]] = true;
  86         }
  87
  88     forAll(boundaryVertex, pI)
  89     {
  90         if( boundaryVertex[pI] && tm.smoothVertex()[pI] )
  91             FatalErrorIn
  92             (
  93                 "void meshOptimizer::untangleMeshFV()"
  94             ) << "Boundary vertex should not be moved!" << abort(FatalError);
  95     }
  96     # endif
  97
  98     label nBadFaces, nGlobalIter(0), nIter;
  99     const label maxNumGlobalIterations(10);
 100
 101     const faceListPMG& faces = mesh_.faces();
 102     boolList changedFace(faces.size(), true);
 103
 104     labelHashSet badFaces;
 105
 106     do
 107     {
 108         nIter = 0;
 109
 110         label minNumBadFaces(10 * faces.size()), minIter(-1);
 111         do
 112         {
 113             nBadFaces = findBadFaces(badFaces, changedFace);
 114
 115             Info << "Iteration " << nIter
 116                 << ". Number of bad faces is " << nBadFaces << endl;
 117
 118             //- perform optimisation
 119             if( nBadFaces == 0 )
 120                 break;
 121
 122             if( nBadFaces < minNumBadFaces )
 123             {
 124                 minNumBadFaces = nBadFaces;
 125                 minIter = nIter;
 126             }
 127
 128             partTetMesh tetMesh(mesh_, badFaces, (nGlobalIter / 5) + 1);
 129
 130             tetMeshOptimisation tmo(tetMesh);
 131
 132             tmo.optimiseUsingKnuppMetric();
 133
 134             tmo.optimiseUsingMeshUntangler();
 135
 136             tmo.optimiseUsingVolumeOptimizer();
 137
 138             tetMesh.updateOrigMesh(&changedFace);
 139
 140         } while( (nIter < minIter+5) && (++nIter < 50) );
 141
 142         if( (nBadFaces == 0) || (++nGlobalIter >= maxNumGlobalIterations) )
 143             break;
 144
 145         // move boundary vertices
 146         nIter = 0;
 147
 148         do
 149         {
 150             nBadFaces = findBadFaces(badFaces, changedFace);
 151
 152             Info << "Iteration " << nIter
 153                 << ". Number of bad faces is " << nBadFaces << endl;
 154
 155             //- perform optimisation
 156             if( nBadFaces == 0 )
 157                 break;
 158
 159             partTetMesh tetMesh(mesh_, badFaces, 0);
 160             tetMeshOptimisation tmo(tetMesh);
 161
 162             if( nGlobalIter < 2 )
 163             {
 164                 //- the point stays in the plane determined by the point normal
 165                 tmo.optimiseBoundaryVolumeOptimizer(true);
 166             }
 167             else if( nGlobalIter < 5 )
 168             {
 169                 //- move points without any constraints on the movement
 170                 tmo.optimiseBoundarySurfaceLaplace();
 171             }
 172             else
 173             {
 174                 //- move boundary points without any constraints
 175                 tmo.optimiseBoundaryVolumeOptimizer(false);
 176             }
 177
 178             tetMesh.updateOrigMesh(&changedFace);
 179
 180         } while( ++nIter < 2 );
 181     }
 182     while( nBadFaces );
 183
 184     Info << "Finished untangling the mesh" << endl;
 185 }
 186
 187 void meshOptimizer::optimizeLowQualityFaces()
 188 {
 189     label nBadFaces, nIter(0);
 190
 191     const faceListPMG& faces = mesh_.faces();
 192     boolList changedFace(faces.size(), true);
 193
 194     label minNumBadFaces(10 * faces.size()), minIter(-1);
 195     do
 196     {
 197         labelHashSet lowQualityFaces;
 198         nBadFaces = findLowQualityFaces(lowQualityFaces, changedFace);
 199
 200         changedFace = false;
 201         forAllConstIter(labelHashSet, lowQualityFaces, it)
 202         changedFace[it.key()] = true;
 203
 204         Info << "Iteration " << nIter
 205             << ". Number of bad faces is " << nBadFaces << endl;
 206
 207         //- perform optimisation
 208         if( nBadFaces == 0 )
 209             break;
 210
 211         if( nBadFaces < minNumBadFaces )
 212         {
 213             minNumBadFaces = nBadFaces;
 214             minIter = nIter;
 215         }
 216
 217         partTetMesh tetMesh(mesh_, lowQualityFaces, 2);
 218
 219         tetMeshOptimisation tmo(tetMesh);
 220
 221         tmo.optimiseUsingKnuppMetric();
 222
 223         tmo.optimiseUsingMeshUntangler();
 224
 225         tmo.optimiseUsingVolumeOptimizer();
 226
 227         tetMesh.updateOrigMesh(&changedFace);
 228
 229     } while( (nIter < minIter+2) && (++nIter < 10) );
 230 }
 231
 232 void meshOptimizer::optimizeMeshFV()
 233 {
 234     Info << "Starting smoothing the mesh" << endl;
 235
 236     laplaceSmoother lps(mesh_, vertexLocation_);
 237     lps.optimizeLaplacianPC(5);
 238
 239     untangleMeshFV();
 240
 241     Info << "Finished smoothing the mesh" << endl;
 242 }
 243
 244 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 245
 246 } // End namespace Foam
 247
 248 // ************************************************************************* //