src/meshLibrary/utilities/smoothers/geometry/meshOptimizer/tetMeshOptimisation/tetMeshOptimisation.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 "tetMeshOptimisation.H"
  30 #include "partTetMesh.H"
  31 #include "meshSurfaceEngine.H"
  32 #include "meshSurfaceEngineModifier.H"
  33
  34 #include "partTetMeshSimplex.H"
  35 #include "meshUntangler.H"
  36 #include "volumeOptimizer.H"
  37 #include "knuppMetric.H"
  38
  39 #include <map>
  40
  41 # ifdef USE_OMP
  42 #include <omp.h>
  43 # endif
  44
  45 // #define DEBUGSearch
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 namespace Foam
  50 {
  51
  52 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  53
  54 // Construct from mesh
  55 tetMeshOptimisation::tetMeshOptimisation(partTetMesh& mesh)
  56 :
  57     tetMesh_(mesh)
  58 {}
  59
  60 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  61
  62 tetMeshOptimisation::~tetMeshOptimisation()
  63 {}
  64
  65 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  66
  67 void tetMeshOptimisation::optimiseUsingKnuppMetric()
  68 {
  69     const LongList<point>& points = tetMesh_.points();
  70     const LongList<partTet>& tets = tetMesh_.tets();
  71     const LongList<direction>& smoothVertex = tetMesh_.smoothVertex();
  72
  73     boolList negativeNode(smoothVertex.size()), invertedTets(tets.size());
  74
  75     //- try getting rid of negative volume using the Patrik Knupp's metric
  76     //- which gets non-negative contributions from invertex tets, only
  77     # ifdef USE_OMP
  78     # pragma omp parallel for if( tets.size() > 100 ) \
  79     schedule(dynamic, 10)
  80     # endif
  81     forAll(tets, tetI)
  82     {
  83         invertedTets[tetI] = false;
  84
  85         if( tets[tetI].mag(points) < VSMALL )
  86             invertedTets[tetI] = true;
  87     }
  88
  89     label nIter(0), nNegative, nNegativeBefore;
  90
  91     do
  92     {
  93         //- find the number of inverted tets
  94         nNegative = 0;
  95         negativeNode = false;
  96         # ifdef USE_OMP
  97         # pragma omp parallel for if( tets.size() > 100 ) \
  98         schedule(dynamic, 10) reduction(+ : nNegative)
  99         # endif
 100         forAll(invertedTets, tetI)
 101         {
 102             if( invertedTets[tetI] )
 103             {
 104                 ++nNegative;
 105                 const partTet& tet = tets[tetI];
 106
 107                 for(label i=0;i<4;++i)
 108                     negativeNode[tet[i]] = true;
 109             }
 110         }
 111
 112         reduce(nNegative, sumOp<label>());
 113         if( nNegative == 0 )
 114             return;
 115
 116         //- make sure that the points at procesor boundaries are selected
 117         //- at all processors
 118         if( Pstream::parRun() )
 119             unifyNegativePoints(negativeNode);
 120
 121         //- smooth the mesh
 122         List<LongList<labelledPoint> > newPositions;
 123         # ifdef USE_OMP
 124         # pragma omp parallel if( smoothVertex.size() > 100 )
 125         # endif
 126         {
 127             # ifdef USE_OMP
 128             # pragma omp master
 129             {
 130                 newPositions.setSize(omp_get_num_threads());
 131             }
 132
 133             # pragma omp barrier
 134
 135             LongList<labelledPoint>& np = newPositions[omp_get_thread_num()];
 136             # else
 137             newPositions.setSize(1);
 138             LongList<labelledPoint>& np = newPositions[0];
 139             # endif
 140
 141             # ifdef USE_OMP
 142             # pragma omp for schedule(dynamic, 10)
 143             # endif
 144             forAll(smoothVertex, nodeI)
 145             {
 146                 if( !negativeNode[nodeI] )
 147                     continue;
 148
 149                 if( smoothVertex[nodeI] & partTetMesh::SMOOTH )
 150                 {
 151                     partTetMeshSimplex simplex(tetMesh_, nodeI);
 152                     knuppMetric(simplex).optimizeNodePosition();
 153                     np.append(labelledPoint(nodeI, simplex.centrePoint()));
 154                 }
 155             }
 156         }
 157
 158         //- update mesh vertices
 159         tetMesh_.updateVerticesSMP(newPositions);
 160         newPositions.clear();
 161
 162         if( Pstream::parRun() )
 163         {
 164             updateBufferLayerPoints();
 165             unifyCoordinatesParallel(&negativeNode);
 166         }
 167
 168         //- check which tets have been repaired
 169         boolList helper(invertedTets.size());
 170         nNegativeBefore = nNegative;
 171         nNegative = 0;
 172
 173         # ifdef USE_OMP
 174         # pragma omp parallel for if( tets.size() > 100 ) \
 175         schedule(dynamic, 10) reduction(+ : nNegative)
 176         # endif
 177         forAll(tets, tetI)
 178         {
 179             helper[tetI] = false;
 180
 181             if( invertedTets[tetI] && (tets[tetI].mag(points) < VSMALL) )
 182             {
 183                 helper[tetI] = true;
 184
 185                 const partTet& tet = tets[tetI];
 186
 187                 for(label i=0;i<4;++i)
 188                     negativeNode[tet[i]] = true;
 189             }
 190         }
 191         invertedTets.transfer(helper);
 192
 193         reduce(nNegative, sumOp<label>());
 194         if( nNegative == 0 )
 195             return;
 196
 197     } while( (nNegative < nNegativeBefore) || (++nIter < 5) );
 198 }
 199
 200 void tetMeshOptimisation::optimiseUsingMeshUntangler()
 201 {
 202     const LongList<point>& points = tetMesh_.points();
 203     const LongList<partTet>& tets = tetMesh_.tets();
 204     const LongList<direction>& smoothVertex = tetMesh_.smoothVertex();
 205
 206     boolList negativeNode(smoothVertex.size()), invertedTets(tets.size());
 207
 208     //- try getting rid of negative volume using the untangler
 209     # ifdef USE_OMP
 210     # pragma omp parallel for if( tets.size() > 100 ) \
 211     schedule(dynamic, 10)
 212     # endif
 213     forAll(tets, tetI)
 214     {
 215         invertedTets[tetI] = false;
 216
 217         if( tets[tetI].mag(points) < VSMALL )
 218             invertedTets[tetI] = true;
 219     }
 220
 221     label nIter(0), nNegative, nNegativeBefore;
 222
 223     do
 224     {
 225         //- find the number of inverted tets
 226         nNegative = 0;
 227         negativeNode = false;
 228         # ifdef USE_OMP
 229         # pragma omp parallel for if( tets.size() > 100 ) \
 230         schedule(dynamic, 10) reduction(+ : nNegative)
 231         # endif
 232         forAll(invertedTets, tetI)
 233         {
 234             if( invertedTets[tetI] )
 235             {
 236                 ++nNegative;
 237                 const partTet& tet = tets[tetI];
 238
 239                 for(label i=0;i<4;++i)
 240                     negativeNode[tet[i]] = true;
 241             }
 242         }
 243
 244         reduce(nNegative, sumOp<label>());
 245         if( nNegative == 0 )
 246             return;
 247
 248         //- make sure that the points at procesor boundaries are selected
 249         //- at all processors
 250         if( Pstream::parRun() )
 251             unifyNegativePoints(negativeNode);
 252
 253         //- smooth the mesh
 254         List<LongList<labelledPoint> > newPositions;
 255         # ifdef USE_OMP
 256         # pragma omp parallel if( smoothVertex.size() > 100 )
 257         # endif
 258         {
 259             # ifdef USE_OMP
 260             # pragma omp master
 261             {
 262                 newPositions.setSize(omp_get_num_threads());
 263             }
 264
 265             # pragma omp barrier
 266
 267             LongList<labelledPoint>& np = newPositions[omp_get_thread_num()];
 268             # else
 269             newPositions.setSize(1);
 270             LongList<labelledPoint>& np = newPositions[0];
 271             # endif
 272
 273             # ifdef USE_OMP
 274             # pragma omp for schedule(dynamic, 10)
 275             # endif
 276             forAll(smoothVertex, nodeI)
 277             {
 278                 if( !negativeNode[nodeI] )
 279                     continue;
 280
 281                 if( smoothVertex[nodeI] & partTetMesh::SMOOTH )
 282                 {
 283                     partTetMeshSimplex simplex(tetMesh_, nodeI);
 284                     meshUntangler(simplex).optimizeNodePosition();
 285                     np.append(labelledPoint(nodeI, simplex.centrePoint()));
 286                 }
 287             }
 288         }
 289
 290         //- update mesh vertices
 291         tetMesh_.updateVerticesSMP(newPositions);
 292         newPositions.clear();
 293
 294         if( Pstream::parRun() )
 295         {
 296             updateBufferLayerPoints();
 297             unifyCoordinatesParallel(&negativeNode);
 298         }
 299
 300         //- check which tets have been repaired
 301         boolList helper(invertedTets.size());
 302         nNegativeBefore = nNegative;
 303         nNegative = 0;
 304         # ifdef USE_OMP
 305         # pragma omp parallel for if( tets.size() > 100 ) \
 306         schedule(dynamic, 10) reduction(+ : nNegative)
 307         # endif
 308         forAll(tets, tetI)
 309         {
 310             helper[tetI] = false;
 311
 312             if( invertedTets[tetI] && (tets[tetI].mag(points) < VSMALL) )
 313             {
 314                 ++nNegative;
 315                 const partTet& tet = tets[tetI];
 316
 317                 for(label i=0;i<4;++i)
 318                     negativeNode[tet[i]] = true;
 319             }
 320         }
 321
 322         reduce(nNegative, sumOp<label>());
 323         if( nNegative == 0 )
 324             return;
 325         invertedTets.transfer(helper);
 326
 327     } while( (nNegative < nNegativeBefore) || (++nIter < 5) );
 328 }
 329
 330 void tetMeshOptimisation::optimiseUsingVolumeOptimizer()
 331 {
 332     const LongList<direction>& smoothVertex = tetMesh_.smoothVertex();
 333
 334     //- use mesh optimizer to improve the result
 335     for(label i=0;i<2;++i)
 336     {
 337         List<LongList<labelledPoint> > newPositions;
 338
 339         # ifdef USE_OMP
 340         # pragma omp parallel if( smoothVertex.size() > 100 )
 341         # endif
 342         {
 343             # ifdef USE_OMP
 344             # pragma omp master
 345             {
 346                 newPositions.setSize(omp_get_num_threads());
 347             }
 348
 349             # pragma omp barrier
 350
 351             LongList<labelledPoint>& np = newPositions[omp_get_thread_num()];
 352             # else
 353             newPositions.setSize(1);
 354             LongList<labelledPoint>& np = newPositions[0];
 355             # endif
 356
 357             # ifdef USE_OMP
 358             # pragma omp for schedule(dynamic, 10)
 359             # endif
 360             forAll(smoothVertex, nodeI)
 361                 if( smoothVertex[nodeI] & partTetMesh::SMOOTH )
 362                 {
 363                     partTetMeshSimplex simplex(tetMesh_, nodeI);
 364
 365                     volumeOptimizer vOpt(simplex);
 366                     vOpt.optimizeNodePosition(1e-5);
 367
 368                     np.append(labelledPoint(nodeI, simplex.centrePoint()));
 369                 }
 370         }
 371
 372         //- update mesh vertices
 373         tetMesh_.updateVerticesSMP(newPositions);
 374         newPositions.clear();
 375
 376         if( Pstream::parRun() )
 377         {
 378             updateBufferLayerPoints();
 379             unifyCoordinatesParallel();
 380         }
 381     }
 382 }
 383
 384 void tetMeshOptimisation::optimiseBoundaryVolumeOptimizer
 385 (
 386     const bool nonShrinking
 387 )
 388 {
 389     const LongList<point>& points = tetMesh_.points();
 390     const LongList<direction>& smoothVertex = tetMesh_.smoothVertex();
 391
 392     # ifdef USE_OMP
 393     label nThreads = omp_get_num_procs();
 394     if( smoothVertex.size() < 100 )
 395       nThreads = 1;
 396     # else
 397     const label nThreads(1);
 398     # endif
 399
 400     for(label i=0;i<3;++i)
 401     {
 402         List<LongList<labelledPoint> > newPositions(nThreads);
 403
 404         # ifdef USE_OMP
 405         # pragma omp parallel num_threads(nThreads)
 406         # endif
 407         {
 408             # ifdef USE_OMP
 409             LongList<labelledPoint>& np = newPositions[omp_get_thread_num()];
 410             # else
 411             LongList<labelledPoint>& np = newPositions[0];
 412             # endif
 413
 414             # ifdef USE_OMP
 415             # pragma omp for schedule(dynamic, 5)
 416             # endif
 417             forAll(smoothVertex, nodeI)
 418                 if( smoothVertex[nodeI] & partTetMesh::BOUNDARY )
 419                 {
 420                     partTetMeshSimplex simplex(tetMesh_, nodeI);
 421
 422                     volumeOptimizer vOpt(simplex);
 423                     vOpt.optimizeNodePosition(1e-5);
 424
 425                     if( nonShrinking )
 426                     {
 427                         //- find boundary faces of the simplex
 428                         const DynList<point, 128>& pts = simplex.pts();
 429                         const DynList<partTet, 128>& tets = simplex.tets();
 430                         DynList<edge, 64> bEdges;
 431                         DynList<label, 64> numAppearances;
 432
 433                         forAll(tets, tetI)
 434                         {
 435                             const partTet& tet = tets[tetI];
 436                             for(label i=0;i<3;++i)
 437                             {
 438                                 edge e(tet[i], tet[(i+1)%3]);
 439
 440                                 const label pos = bEdges.containsAtPosition(e);
 441                                 if( pos < 0 )
 442                                 {
 443                                     bEdges.append(e);
 444                                     numAppearances.append(1);
 445                                 }
 446                                 else
 447                                 {
 448                                     ++numAppearances(pos);
 449                                 }
 450                             }
 451                         }
 452
 453                         //- create normal tensor of the simplex
 454                         symmTensor nt(symmTensor::zero);
 455                         forAll(bEdges, beI)
 456                         {
 457                             if( numAppearances[beI] != 1 )
 458                                 continue;
 459
 460                             triangle<point, point> tri
 461                             (
 462                                 pts[bEdges[beI].start()],
 463                                 pts[bEdges[beI].end()],
 464                                 points[nodeI]
 465                             );
 466
 467                             vector n = tri.normal();
 468                             n /= (mag(n) + VSMALL);
 469                             for(direction i=0;i<vector::nComponents;++i)
 470                                 if( Foam::mag(n[i]) < (100. * SMALL) )
 471                                     n[i] = 0.0;
 472
 473                             nt += symm(n * n);
 474                         }
 475
 476                         const vector ev = eigenValues(nt);
 477
 478                         //- make sure the point stays on the surface
 479                         vector disp = simplex.centrePoint() - points[nodeI];
 480
 481                         if( mag(ev[2]) > (mag(ev[1]) + mag(ev[0])) )
 482                         {
 483                             //- ordinary surface vertex
 484                             vector normal = eigenVector(nt, ev[2]);
 485                             normal /= (mag(normal)+VSMALL);
 486                             disp -= (disp & normal) * normal;
 487                         }
 488                         else if( mag(ev[1]) > 0.5 * (mag(ev[2]) + mag(ev[0])) )
 489                         {
 490                             //- this vertex is on an edge
 491                             vector normal1 = eigenVector(nt, ev[1]);
 492                             normal1 /= (mag(normal1)+VSMALL);
 493                             vector normal2 = eigenVector(nt, ev[2]);
 494                             normal2 /= (mag(normal2)+VSMALL);
 495
 496                             vector eVec = normal1 ^ normal2;
 497                             eVec /= (mag(eVec) + VSMALL);
 498
 499                             disp = (disp & eVec) * eVec;
 500                         }
 501                         else
 502                         {
 503                             //- this vertex is a corner. do not move it
 504                             continue;
 505                         }
 506
 507                         const point newP = points[nodeI] + disp;
 508                         np.append(labelledPoint(nodeI, newP));
 509                     }
 510                     else
 511                     {
 512                         //- move the vertex without constraining it
 513                         np.append(labelledPoint(nodeI, simplex.centrePoint()));
 514                     }
 515                 }
 516         }
 517
 518         //- update tetMesh
 519         tetMesh_.updateVerticesSMP(newPositions);
 520         newPositions.clear();
 521
 522         if( Pstream::parRun() )
 523         {
 524             updateBufferLayerPoints();
 525             unifyCoordinatesParallel();
 526         }
 527     }
 528 }
 529
 530 void tetMeshOptimisation::optimiseBoundarySurfaceLaplace()
 531 {
 532     const LongList<direction>& smoothVertex = tetMesh_.smoothVertex();
 533
 534     # ifdef USE_OMP
 535     label nThreads = omp_get_num_procs();
 536     if( smoothVertex.size() < 1000 )
 537       nThreads = 1;
 538     # else
 539     const label nThreads(1);
 540     # endif
 541
 542     for(label i=0;i<3;++i)
 543     {
 544         List<LongList<labelledPoint> > newPositions(nThreads);
 545
 546         # ifdef USE_OMP
 547         # pragma omp parallel num_threads(nThreads)
 548         # endif
 549         {
 550             # ifdef USE_OMP
 551             LongList<labelledPoint>& np =
 552                 newPositions[omp_get_thread_num()];
 553             # else
 554             LongList<labelledPoint>& np = newPositions[0];
 555             # endif
 556
 557             # ifdef USE_OMP
 558             # pragma omp for schedule(dynamic, 5)
 559             # endif
 560             forAll(smoothVertex, nodeI)
 561             {
 562                 if( smoothVertex[nodeI] & partTetMesh::BOUNDARY )
 563                 {
 564                     partTetMeshSimplex simplex(tetMesh_, nodeI);
 565
 566                     //- find boundary faces of the simplex
 567                     const DynList<point, 128>& pts = simplex.pts();
 568                     const DynList<partTet, 128>& tets = simplex.tets();
 569                     DynList<edge, 64> bndEdges;
 570                     DynList<label, 64> numAppearances;
 571
 572                     //- find boundary edges of the simplex
 573                     forAll(tets, tetI)
 574                     {
 575                         const partTet& tet = tets[tetI];
 576                         for(label i=0;i<3;++i)
 577                         {
 578                             const edge e(tet[i], tet[(i+1)%3]);
 579                             const label pos = bndEdges.containsAtPosition(e);
 580
 581                             if( pos < 0 )
 582                             {
 583                                 bndEdges.append(e);
 584                                 numAppearances.append(1);
 585                             }
 586                             else
 587                             {
 588                                 ++numAppearances(pos);
 589                             }
 590                         }
 591                     }
 592
 593                     point newP(vector::zero);
 594                     label counter(0);
 595                     forAll(bndEdges, beI)
 596                     {
 597                         if( numAppearances[beI] != 1 )
 598                             continue;
 599
 600                         triangle<point, point> tri
 601                         (
 602                             pts[bndEdges[beI].start()],
 603                             pts[bndEdges[beI].end()],
 604                             simplex.centrePoint()
 605                         );
 606
 607                         newP += tri.centre();
 608                         ++counter;
 609                     }
 610
 611                     if( counter != 0 )
 612                     {
 613                         newP /= counter;
 614                         np.append(labelledPoint(nodeI, newP));
 615                     }
 616                 }
 617             }
 618         }
 619
 620         //- update tetMesh with new vertex positions
 621         tetMesh_.updateVerticesSMP(newPositions);
 622         newPositions.clear();
 623
 624         if( Pstream::parRun() )
 625         {
 626             updateBufferLayerPoints();
 627             unifyCoordinatesParallel();
 628         }
 629     }
 630 }
 631
 632 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 633
 634 } // End namespace Foam
 635
 636 // ************************************************************************* //