optimizeMeshFV.C 18.2 KB
Newer Older
franjo's avatar
franjo committed
1
2
/*---------------------------------------------------------------------------*\
  =========                 |
Tomislav Lugaric's avatar
Tomislav Lugaric committed
3
  \\      /  F ield         | cfMesh: A library for mesh generation
franjo's avatar
franjo committed
4
   \\    /   O peration     |
Tomislav Lugaric's avatar
Tomislav Lugaric committed
5
6
    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
franjo's avatar
franjo committed
7
8
-------------------------------------------------------------------------------
License
Tomislav Lugaric's avatar
Tomislav Lugaric committed
9
    This file is part of cfMesh.
franjo's avatar
franjo committed
10

Tomislav Lugaric's avatar
Tomislav Lugaric committed
11
    cfMesh is free software; you can redistribute it and/or modify it
franjo's avatar
franjo committed
12
    under the terms of the GNU General Public License as published by the
Tomislav Lugaric's avatar
Tomislav Lugaric committed
13
    Free Software Foundation; either version 3 of the License, or (at your
franjo's avatar
franjo committed
14
15
    option) any later version.

Tomislav Lugaric's avatar
Tomislav Lugaric committed
16
    cfMesh is distributed in the hope that it will be useful, but WITHOUT
franjo's avatar
franjo committed
17
18
19
20
21
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
Tomislav Lugaric's avatar
Tomislav Lugaric committed
22
    along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
franjo's avatar
franjo committed
23
24
25
26
27

Description

\*---------------------------------------------------------------------------*/

28
29
#include <stdexcept>

franjo's avatar
franjo committed
30
31
32
#include "demandDrivenData.H"
#include "meshOptimizer.H"
#include "polyMeshGenAddressing.H"
33
#include "polyMeshGenChecks.H"
franjo's avatar
franjo committed
34
35
36
37
#include "partTetMesh.H"
#include "HashSet.H"

#include "tetMeshOptimisation.H"
38
#include "boundaryLayerOptimisation.H"
Franjo's avatar
Franjo committed
39
#include "refineBoundaryLayers.H"
40
#include "meshSurfaceEngine.H"
franjo's avatar
franjo committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

//#define DEBUGSmooth

# ifdef DEBUGSmooth
#include "helperFunctions.H"
#include "polyMeshGenModifier.H"
# endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Franjo's avatar
Franjo committed
56
57
58
59
void meshOptimizer::untangleMeshFV
(
    const label maxNumGlobalIterations,
    const label maxNumIterations,
60
61
    const label maxNumSurfaceIterations,
    const bool relaxedCheck
Franjo's avatar
Franjo committed
62
)
franjo's avatar
franjo committed
63
{
64
    Info << "Starting untangling the mesh" << endl;
65

66
67
68
69
70
71
72
    # ifdef DEBUGSmooth
    partTetMesh tm(mesh_);
    forAll(tm.tets(), tetI)
        if( tm.tets()[tetI].mag(tm.points()) < 0.0 )
            Info << "Tet " << tetI << " is inverted!" << endl;
    polyMeshGen tetPolyMesh(mesh_.returnTime());
    tm.createPolyMesh(tetPolyMesh);
franjo's avatar
franjo committed
73
    polyMeshGenModifier(tetPolyMesh).removeUnusedVertices();
74
75
76
    forAll(tm.smoothVertex(), pI)
        if( !tm.smoothVertex()[pI] )
            Info << "Point " << pI << " cannot be moved!" << endl;
77

franjo's avatar
franjo committed
78
79
80
81
82
83
84
    const VRWGraph& pTets = tm.pointTets();
    forAll(pTets, pointI)
    {
        const LongList<partTet>& tets = tm.tets();
        forAllRow(pTets, pointI, i)
            if( tets[pTets(pointI, i)].whichPosition(pointI) < 0 )
                FatalError << "Wrong partTet" << abort(FatalError);
85

franjo's avatar
franjo committed
86
87
        partTetMeshSimplex simplex(tm, pointI);
    }
88

89
90
91
92
93
94
    boolList boundaryVertex(tetPolyMesh.points().size(), false);
    const labelList& neighbour = tetPolyMesh.neighbour();
    forAll(neighbour, faceI)
        if( neighbour[faceI] == -1 )
        {
            const face& f = tetPolyMesh.faces()[faceI];
95

96
97
98
            forAll(f, pI)
                boundaryVertex[f[pI]] = true;
        }
99

100
101
102
103
104
105
106
107
108
    forAll(boundaryVertex, pI)
    {
        if( boundaryVertex[pI] && tm.smoothVertex()[pI] )
            FatalErrorIn
            (
                "void meshOptimizer::untangleMeshFV()"
            ) << "Boundary vertex should not be moved!" << abort(FatalError);
    }
    # endif
109

110
    label nBadFaces, nGlobalIter(0), nIter;
111

112
    const faceListPMG& faces = mesh_.faces();
113

114
    boolList changedFace(faces.size(), true);
115

Franjo's avatar
Franjo committed
116
117
118
119
120
121
122
    //- check if any points in the tet mesh shall not move
    labelLongList lockedPoints;
    forAll(vertexLocation_, pointI)
    {
        if( vertexLocation_[pointI] & LOCKED )
            lockedPoints.append(pointI);
    }
123

franjo's avatar
franjo committed
124
    labelHashSet badFaces;
125

126
127
128
    do
    {
        nIter = 0;
129

franjo's avatar
franjo committed
130
        label minNumBadFaces(10 * faces.size()), minIter(-1);
131
132
        do
        {
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
            if( !relaxedCheck )
            {
                nBadFaces =
                    polyMeshGenChecks::findBadFaces
                    (
                        mesh_,
                        badFaces,
                        false,
                        &changedFace
                    );
            }
            else
            {
                nBadFaces =
                    polyMeshGenChecks::findBadFacesRelaxed
                    (
                        mesh_,
                        badFaces,
                        false,
                        &changedFace
                    );
            }
155

franjo's avatar
franjo committed
156
157
            Info << "Iteration " << nIter
                << ". Number of bad faces is " << nBadFaces << endl;
158

159
160
161
            //- perform optimisation
            if( nBadFaces == 0 )
                break;
162

franjo's avatar
franjo committed
163
164
165
166
167
            if( nBadFaces < minNumBadFaces )
            {
                minNumBadFaces = nBadFaces;
                minIter = nIter;
            }
168

169
            //- create a tet mesh from the mesh and the labels of bad faces
Franjo's avatar
Franjo committed
170
171
172
173
174
175
176
            partTetMesh tetMesh
            (
                mesh_,
                lockedPoints,
                badFaces,
                (nGlobalIter / 2) + 1
            );
franjo's avatar
franjo committed
177

178
179
            //- construct tetMeshOptimisation and improve positions of
            //- points in the tet mesh
franjo's avatar
franjo committed
180
            tetMeshOptimisation tmo(tetMesh);
181

franjo's avatar
franjo committed
182
            tmo.optimiseUsingKnuppMetric();
183

franjo's avatar
franjo committed
184
            tmo.optimiseUsingMeshUntangler();
185

franjo's avatar
franjo committed
186
            tmo.optimiseUsingVolumeOptimizer();
187

188
            //- update points in the mesh from the coordinates in the tet mesh
189
            tetMesh.updateOrigMesh(&changedFace);
190

Franjo's avatar
Franjo committed
191
        } while( (nIter < minIter+5) && (++nIter < maxNumIterations) );
franjo's avatar
franjo committed
192

193
194
        if( (nBadFaces == 0) || (++nGlobalIter >= maxNumGlobalIterations) )
            break;
195

franjo's avatar
franjo committed
196
197
        // move boundary vertices
        nIter = 0;
198

Franjo's avatar
Franjo committed
199
        while( nIter++ < maxNumSurfaceIterations );
200
        {
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
            if( !relaxedCheck )
            {
                nBadFaces =
                    polyMeshGenChecks::findBadFaces
                    (
                        mesh_,
                        badFaces,
                        false,
                        &changedFace
                    );
            }
            else
            {
                nBadFaces =
                    polyMeshGenChecks::findBadFacesRelaxed
                    (
                        mesh_,
                        badFaces,
                        false,
                        &changedFace
                    );
            }
223

franjo's avatar
franjo committed
224
225
            Info << "Iteration " << nIter
                << ". Number of bad faces is " << nBadFaces << endl;
226

227
228
            //- perform optimisation
            if( nBadFaces == 0 )
229
            {
230
                break;
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
            }
            else if( enforceConstraints_ )
            {
                const label subsetId =
                    mesh_.addPointSubset(badPointsSubsetName_);

                forAllConstIter(labelHashSet, badFaces, it)
                {
                    const face& f = faces[it.key()];
                    forAll(f, pI)
                        mesh_.addPointToSubset(subsetId, f[pI]);
                }

                WarningIn
                (
                    "void meshOptimizer::untangleMeshFV()"
                ) << "Writing mesh with " << badPointsSubsetName_
                  << " subset. These points cannot be untangled"
                  << " without sacrificing geometry constraints. Exitting.."
                  << endl;

                returnReduce(1, sumOp<label>());
                throw std::logic_error
                (
                    "void meshOptimizer::untangleMeshFV()"
                    "Cannot untangle mesh!!"
                );
            }
259

Franjo's avatar
Franjo committed
260
261
            //- create tethrahedral mesh from the cells which shall be smoothed
            partTetMesh tetMesh(mesh_, lockedPoints, badFaces, 0);
262
263

            //- contruct tetMeshOptimisation
franjo's avatar
franjo committed
264
            tetMeshOptimisation tmo(tetMesh);
265

franjo's avatar
franjo committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
            if( nGlobalIter < 2 )
            {
                //- the point stays in the plane determined by the point normal
                tmo.optimiseBoundaryVolumeOptimizer(true);
            }
            else if( nGlobalIter < 5 )
            {
                //- move points without any constraints on the movement
                tmo.optimiseBoundarySurfaceLaplace();
            }
            else
            {
                //- move boundary points without any constraints
                tmo.optimiseBoundaryVolumeOptimizer(false);
            }
281

282
            tetMesh.updateOrigMesh(&changedFace);
283

Franjo's avatar
Franjo committed
284
        }
285
286

    } while( nBadFaces );
287

288
289
290
291
292
293
294
295
296
297
298
    if( nBadFaces != 0 )
    {
        label subsetId = mesh_.faceSubsetIndex("badFaces");
        if( subsetId >= 0 )
            mesh_.removeFaceSubset(subsetId);
        subsetId = mesh_.addFaceSubset("badFaces");

        forAllConstIter(labelHashSet, badFaces, it)
            mesh_.addFaceToSubset(subsetId, it.key());
    }

299
    Info << "Finished untangling the mesh" << endl;
franjo's avatar
franjo committed
300
301
}

302
void meshOptimizer::optimizeBoundaryLayer(const bool addBufferLayer)
303
{
304
305
306
307
308
    if( mesh_.returnTime().foundObject<IOdictionary>("meshDict") )
    {
        const dictionary& meshDict =
            mesh_.returnTime().lookupObject<IOdictionary>("meshDict");

309
310
        bool smoothLayer(false);

311
312
313
314
315
        if( meshDict.found("boundaryLayers") )
        {
            const dictionary& layersDict = meshDict.subDict("boundaryLayers");

            if( layersDict.found("optimiseLayer") )
316
                smoothLayer = readBool(layersDict.lookup("optimiseLayer"));
317
        }
318

319
320
321
        if( !smoothLayer )
            return;

322
        if( addBufferLayer )
Franjo's avatar
Franjo committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
        {
            //- create a buffer layer which will not be modified by the smoother
            refineBoundaryLayers refLayers(mesh_);

            refineBoundaryLayers::readSettings(meshDict, refLayers);

            refLayers.activateSpecialMode();

            refLayers.refineLayers();

            clearSurface();
            calculatePointLocations();
        }

337
        Info << "Starting optimising boundary layer" << endl;
338

339
340
        const meshSurfaceEngine& mse = meshSurface();
        const labelList& faceOwner = mse.faceOwners();
341

342
        boundaryLayerOptimisation optimiser(mesh_, mse);
Franjo's avatar
Franjo committed
343

344
        boundaryLayerOptimisation::readSettings(meshDict, optimiser);
345

346
        optimiser.optimiseLayer();
347

348
        //- check if the bnd layer is tangled somewhere
349
        labelLongList bndLayerCells;
350
351
        const boolList& baseFace = optimiser.isBaseFace();

352
353
354
355
        # ifdef DEBUGSmooth
        const label blCellsId = mesh_.addCellSubset("blCells");
        # endif

356
        forAll(baseFace, bfI)
Franjo's avatar
Franjo committed
357
        {
358
            if( baseFace[bfI] )
359
            {
360
                bndLayerCells.append(faceOwner[bfI]);
361
362
363
364
365

                # ifdef DEBUGSmooth
                mesh_.addCellToSubset(blCellsId, faceOwner[bfI]);
                # endif
            }
Franjo's avatar
Franjo committed
366
        }
367

368
369
        clearSurface();
        mesh_.clearAddressingData();
370

371
372
        //- lock boundary layer points, faces and cells
        lockCells(bndLayerCells);
Franjo's avatar
Franjo committed
373

374
375
376
377
378
379
380
        # ifdef DEBUGSmooth
        pointField origPoints(mesh_.points().size());
        forAll(origPoints, pI)
            origPoints[pI] = mesh_.points()[pI];
        # endif

        //- optimize mesh quality
Franjo's avatar
Franjo committed
381
        optimizeMeshFV(5, 1, 50, 0);
382

383
384
        //- untangle remaining faces and lock the boundary layer cells
        untangleMeshFV(2, 50, 0);
385

386
387
388
389
390
391
392
393
394
395
396
397
        # ifdef DEBUGSmooth
        forAll(vertexLocation_, pI)
        {
            if( vertexLocation_[pI] & LOCKED )
            {
                if( mag(origPoints[pI] - mesh_.points()[pI]) > SMALL )
                    FatalError << "Locked points were moved"
                               << abort(FatalError);
            }
        }
        # endif

398
399
400
401
402
        //- unlock bnd layer points
        removeUserConstraints();

        Info << "Finished optimising boundary layer" << endl;
    }
403
404
}

405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
void meshOptimizer::untangleBoundaryLayer()
{
    bool untangleLayer(true);
    if( mesh_.returnTime().foundObject<IOdictionary>("meshDict") )
    {
        const dictionary& meshDict =
            mesh_.returnTime().lookupObject<IOdictionary>("meshDict");

        if( meshDict.found("boundaryLayers") )
        {
            const dictionary& layersDict = meshDict.subDict("boundaryLayers");

            if( layersDict.found("untangleLayers") )
            {
                untangleLayer =
                    readBool(layersDict.lookup("untangleLayers"));
            }
        }
    }

    if( !untangleLayer )
    {
        labelHashSet badFaces;
428
        polyMeshGenChecks::checkFacePyramids(mesh_, false, VSMALL, &badFaces);
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451

        const label nInvalidFaces =
            returnReduce(badFaces.size(), sumOp<label>());

        if( nInvalidFaces != 0 )
        {
            const labelList& owner = mesh_.owner();
            const labelList& neighbour = mesh_.neighbour();

            const label badBlCellsId =
                mesh_.addCellSubset("invalidBoundaryLayerCells");

            forAllConstIter(labelHashSet, badFaces, it)
            {
                mesh_.addCellToSubset(badBlCellsId, owner[it.key()]);

                if( neighbour[it.key()] < 0 )
                    continue;

                mesh_.addCellToSubset(badBlCellsId, neighbour[it.key()]);
            }

            returnReduce(1, sumOp<label>());
452
453

            throw std::logic_error
454
455
            (
                "void meshOptimizer::untangleBoundaryLayer()"
456
457
458
                "Found invalid faces in the boundary layer."
                " Cannot untangle mesh!!"
            );
459
460
461
462
        }
    }
    else
    {
463
        optimizeLowQualityFaces();
464
        removeUserConstraints();
465
        untangleMeshFV(2, 50, 1, true);
466
467
468
    }
}

Franjo's avatar
Franjo committed
469
void meshOptimizer::optimizeLowQualityFaces(const label maxNumIterations)
franjo's avatar
franjo committed
470
471
{
    label nBadFaces, nIter(0);
472

473
474
    const faceListPMG& faces = mesh_.faces();
    boolList changedFace(faces.size(), true);
475

Franjo's avatar
Franjo committed
476
477
478
479
480
481
482
    //- check if any points in the tet mesh shall not move
    labelLongList lockedPoints;
    forAll(vertexLocation_, pointI)
    {
        if( vertexLocation_[pointI] & LOCKED )
            lockedPoints.append(pointI);
    }
483

franjo's avatar
franjo committed
484
485
486
    do
    {
        labelHashSet lowQualityFaces;
487
488
489
490
491
492
493
494
        nBadFaces =
            polyMeshGenChecks::findLowQualityFaces
            (
                mesh_,
                lowQualityFaces,
                false,
                &changedFace
            );
495

franjo's avatar
franjo committed
496
497
        changedFace = false;
        forAllConstIter(labelHashSet, lowQualityFaces, it)
498
            changedFace[it.key()] = true;
499

franjo's avatar
franjo committed
500
501
        Info << "Iteration " << nIter
            << ". Number of bad faces is " << nBadFaces << endl;
502

franjo's avatar
franjo committed
503
504
505
        //- perform optimisation
        if( nBadFaces == 0 )
            break;
506

Franjo's avatar
Franjo committed
507
        partTetMesh tetMesh(mesh_, lockedPoints, lowQualityFaces, 2);
508
509
510

        //- construct tetMeshOptimisation and improve positions
        //- of points in the tet mesh
franjo's avatar
franjo committed
511
        tetMeshOptimisation tmo(tetMesh);
512

franjo's avatar
franjo committed
513
514
        tmo.optimiseUsingVolumeOptimizer();

515
        //- update points in the mesh from the new coordinates in the tet mesh
franjo's avatar
franjo committed
516
517
        tetMesh.updateOrigMesh(&changedFace);

Franjo's avatar
Franjo committed
518
    } while( ++nIter < maxNumIterations );
Franjo's avatar
Franjo committed
519
520
521
522
523
524
525
526
527
528
529
530
531
}

void meshOptimizer::optimizeMeshNearBoundaries
(
    const label maxNumIterations,
    const label numLayersOfCells
)
{
    label nIter(0);

    const faceListPMG& faces = mesh_.faces();
    boolList changedFace(faces.size(), true);

Franjo's avatar
Franjo committed
532
533
534
535
536
537
538
539
540
    //- check if any points in the tet mesh shall not move
    labelLongList lockedPoints;
    forAll(vertexLocation_, pointI)
    {
        if( vertexLocation_[pointI] & LOCKED )
            lockedPoints.append(pointI);
    }

    partTetMesh tetMesh(mesh_, lockedPoints, numLayersOfCells);
Franjo's avatar
Franjo committed
541
542
543
544
    tetMeshOptimisation tmo(tetMesh);
    Info << "Iteration:" << flush;
    do
    {
545
        tmo.optimiseUsingVolumeOptimizer(1);
Franjo's avatar
Franjo committed
546
547
548
549
550
551
552
553

        tetMesh.updateOrigMesh(&changedFace);

        Info << "." << flush;

    } while( ++nIter < maxNumIterations );

    Info << endl;
franjo's avatar
franjo committed
554
}
555

556
557
558
559
560
561
562
void meshOptimizer::optimizeMeshFV
(
    const label numLaplaceIterations,
    const label maxNumGlobalIterations,
    const label maxNumIterations,
    const label maxNumSurfaceIterations
)
franjo's avatar
franjo committed
563
564
{
    Info << "Starting smoothing the mesh" << endl;
565

franjo's avatar
franjo committed
566
    laplaceSmoother lps(mesh_, vertexLocation_);
567
568
569
570
571
572
573
574
    lps.optimizeLaplacianPC(numLaplaceIterations);

    untangleMeshFV
    (
        maxNumGlobalIterations,
        maxNumIterations,
        maxNumSurfaceIterations
    );
franjo's avatar
franjo committed
575
576
577
578

    Info << "Finished smoothing the mesh" << endl;
}

579
580
void meshOptimizer::optimizeMeshFVBestQuality
(
Franjo's avatar
Franjo committed
581
582
    const label maxNumIterations,
    const scalar threshold
583
584
585
)
{
    label nBadFaces, nIter(0);
Franjo's avatar
Cleanup    
Franjo committed
586
    label minIter(-1);
587
588
589
590
591
592
593
594
595
596
597
598

    const faceListPMG& faces = mesh_.faces();
    boolList changedFace(faces.size(), true);

    //- check if any points in the tet mesh shall not move
    labelLongList lockedPoints;
    forAll(vertexLocation_, pointI)
    {
        if( vertexLocation_[pointI] & LOCKED )
            lockedPoints.append(pointI);
    }

Franjo's avatar
Franjo committed
599
    label minNumBadFaces(10 * faces.size());
600
601
602
603
604
605
606
607
608
609
    do
    {
        labelHashSet lowQualityFaces;
        nBadFaces =
            polyMeshGenChecks::findWorstQualityFaces
            (
                mesh_,
                lowQualityFaces,
                false,
                &changedFace,
Franjo's avatar
Franjo committed
610
                threshold
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
            );

        changedFace = false;
        forAllConstIter(labelHashSet, lowQualityFaces, it)
            changedFace[it.key()] = true;

        Info << "Iteration " << nIter
            << ". Number of worst quality faces is " << nBadFaces << endl;

        //- perform optimisation
        if( nBadFaces == 0 )
            break;

        if( nBadFaces < minNumBadFaces )
        {
            minNumBadFaces = nBadFaces;
Franjo's avatar
Cleanup    
Franjo committed
627
628

            //- update the iteration number when the minimum is achieved
629
630
631
632
633
634
635
636
637
            minIter = nIter;
        }

        partTetMesh tetMesh(mesh_, lockedPoints, lowQualityFaces, 2);

        //- construct tetMeshOptimisation and improve positions
        //- of points in the tet mesh
        tetMeshOptimisation tmo(tetMesh);

Franjo's avatar
Franjo committed
638
        tmo.optimiseUsingVolumeOptimizer(20);
639
640
641
642

        //- update points in the mesh from the new coordinates in the tet mesh
        tetMesh.updateOrigMesh(&changedFace);

643
    } while( (nIter < minIter+5) && (++nIter < maxNumIterations) );
644
645
}

franjo's avatar
franjo committed
646
647
648
649
650
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //