src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/boundaryLayerOptimisation/boundaryLayerOptimisationFunctions.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 "boundaryLayerOptimisation.H"
  30 #include "meshSurfacePartitioner.H"
  31 #include "meshSurfaceEngine.H"
  32 #include "detectBoundaryLayers.H"
  33 //#include "helperFunctions.H"
  34 //#include "labelledScalar.H"
  35 //#include "refLabelledPoint.H"
  36 //#include "refLabelledPointScalar.H"
  37 #include "polyMeshGenAddressing.H"
  38 #include "meshSurfaceOptimizer.H"
  39 #include "meshSurfaceEngineModifier.H"
  40 //#include "partTetMeshSimplex.H"
  41 //#include "volumeOptimizer.H"
  42 #include "OFstream.H"
  43
  44 //#define DEBUGLayer
  45
  46 # ifdef USE_OMP
  47 #include <omp.h>
  48 # endif
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 namespace Foam
  53 {
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57 void boundaryLayerOptimisation::writeVTK
  58 (
  59     const fileName& fName,
  60     const pointField& origin,
  61     const vectorField& vecs
  62 )
  63 {
  64     if( origin.size() != vecs.size() )
  65         FatalErrorIn
  66         (
  67             "void boundaryLayerOptimisation::writeVTK(const fileName&,"
  68             " const pointField&, const vectorField&)"
  69         ) << "Sizes do not match" << abort(FatalError);
  70
  71     OFstream file(fName);
  72
  73     //- write the header
  74     file << "# vtk DataFile Version 3.0\n";
  75     file << "vtk output\n";
  76     file << "ASCII\n";
  77     file << "DATASET POLYDATA\n";
  78
  79     //- write points
  80     file << "POINTS " << 2*origin.size() << " float\n";
  81     forAll(origin, pI)
  82     {
  83         const point& p = origin[pI];
  84
  85         file << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
  86
  87         const point op = p + vecs[pI];
  88
  89         file << op.x() << ' ' << op.y() << ' ' << op.z() << nl;
  90     }
  91
  92     //- write lines
  93     file << "\nLINES " << vecs.size()
  94          << " " << 3*vecs.size() << nl;
  95     forAll(vecs, eI)
  96     {
  97         file << 2 << " " << 2*eI << " " << (2*eI+1) << nl;
  98     }
  99
 100     file << "\n";
 101 }
 102
 103 void boundaryLayerOptimisation::writeHairEdges
 104 (
 105     const fileName& fName,
 106     const direction eType,
 107     const vectorField& vecs
 108 ) const
 109 {
 110     if( vecs.size() != hairEdges_.size() )
 111         FatalErrorIn
 112         (
 113             "void boundaryLayerOptimisation::writeHairEdges"
 114             "(const fileName&, const direction, const vectorField&) const"
 115         ) << "Sizes do not match" << abort(FatalError);
 116
 117     //- count the number of hair edges matching this criteria
 118     label counter(0);
 119
 120     forAll(hairEdgeType_, heI)
 121         if( hairEdgeType_[heI] & eType )
 122             ++counter;
 123
 124     //- copy edge vector
 125     vectorField copyVecs(counter);
 126     pointField pts(counter);
 127
 128     counter = 0;
 129
 130     const pointFieldPMG& points = mesh_.points();
 131
 132     forAll(hairEdgeType_, heI)
 133     {
 134         if( hairEdgeType_[heI] & eType )
 135         {
 136             const edge& he = hairEdges_[heI];
 137
 138             pts[counter] = points[he.start()];
 139             copyVecs[counter] = vecs[heI] * he.mag(points);
 140
 141             ++counter;
 142         }
 143     }
 144
 145     //- write data to file
 146     writeVTK(fName, pts, copyVecs);
 147 }
 148
 149 void boundaryLayerOptimisation::writeHairEdges
 150 (
 151     const fileName& fName,
 152     const direction eType
 153 ) const
 154 {
 155     //- count the number of hair edges matching this criteria
 156     label counter(0);
 157
 158     forAll(hairEdgeType_, heI)
 159         if( hairEdgeType_[heI] & eType )
 160             ++counter;
 161
 162     //- copy edge vector
 163     vectorField vecs(counter);
 164     pointField pts(counter);
 165
 166     counter = 0;
 167
 168     const pointFieldPMG& points = mesh_.points();
 169
 170     forAll(hairEdgeType_, heI)
 171     {
 172         if( hairEdgeType_[heI] & eType )
 173         {
 174             const edge& he = hairEdges_[heI];
 175
 176             pts[counter] = points[he.start()];
 177             vecs[counter] = he.vec(points);
 178
 179             ++counter;
 180         }
 181     }
 182
 183     //- write data to file
 184     writeVTK(fName,pts, vecs);
 185 }
 186
 187 const meshSurfaceEngine& boundaryLayerOptimisation::meshSurface() const
 188 {
 189     if( !meshSurfacePtr_ )
 190     {
 191         # ifdef USE_OMP
 192         if( omp_in_parallel() )
 193             FatalErrorIn
 194             (
 195                 "const meshSurfaceEngine&"
 196                 " boundaryLayerOptimisation::meshSurface()"
 197             ) << "Cannot generate meshSurfaceEngine" << abort(FatalError);
 198         # endif
 199
 200         meshSurfacePtr_ = new meshSurfaceEngine(mesh_);
 201     }
 202
 203     return *meshSurfacePtr_;
 204 }
 205
 206 const meshSurfacePartitioner&
 207 boundaryLayerOptimisation::surfacePartitioner() const
 208 {
 209     if( !partitionerPtr_ )
 210     {
 211         # ifdef USE_OMP
 212         if( omp_in_parallel() )
 213             FatalErrorIn
 214             (
 215                 "const meshSurfacePartitioner& "
 216                 "boundaryLayerOptimisation::surfacePartitioner()"
 217             ) << "Cannot generate meshSurfacePartitioner" << abort(FatalError);
 218         # endif
 219
 220         partitionerPtr_ = new meshSurfacePartitioner(meshSurface());
 221     }
 222
 223     return *partitionerPtr_;
 224 }
 225
 226 void boundaryLayerOptimisation::calculateHairEdges()
 227 {
 228     const meshSurfaceEngine& mse = meshSurface();
 229     const edgeList& edges = mse.edges();
 230     const VRWGraph& edgeFaces = mse.edgeFaces();
 231     const VRWGraph& bpEdges = mse.boundaryPointEdges();
 232     const labelList& faceOwner = mse.faceOwners();
 233     const faceList::subList& bFaces = mse.boundaryFaces();
 234     const labelList& bp = mse.bp();
 235
 236     const meshSurfacePartitioner& mPart = surfacePartitioner();
 237
 238     //- detect layers in the mesh
 239     const detectBoundaryLayers detectLayers(mPart);
 240
 241     hairEdges_ = detectLayers.hairEdges();
 242     hairEdgesAtBndPoint_ = detectLayers.hairEdgesAtBndPoint();
 243
 244     //- mark boundary faces which are base face for the boundary layer
 245     const labelList& layerAtBndFace = detectLayers.faceInLayer();
 246     isBndLayerBase_.setSize(bFaces.size());
 247     # ifdef USE_OMP
 248     # pragma omp parallel for schedule(dynamic, 100)
 249     # endif
 250     forAll(layerAtBndFace, bfI)
 251     {
 252         if( layerAtBndFace[bfI] < 0 )
 253         {
 254             isBndLayerBase_[bfI] = false;
 255         }
 256         else
 257         {
 258             isBndLayerBase_[bfI] = true;
 259         }
 260     }
 261
 262     # ifdef DEBUGLayer
 263     const label bndLayerFaceId = mesh_.addFaceSubset("bndLayerFaces");
 264     const label startBndFaceI = mesh_.boundaries()[0].patchStart();
 265     forAll(isBndLayerBase_, bfI)
 266         if( isBndLayerBase_[bfI] )
 267             mesh_.addFaceToSubset(bndLayerFaceId, startBndFaceI+bfI);
 268     # endif
 269
 270     //- check if a face is an exiting face for a bnd layer
 271     isExitFace_.setSize(isBndLayerBase_.size());
 272     isExitFace_ = false;
 273
 274     # ifdef USE_OMP
 275     # pragma omp parallel for schedule(dynamic, 100)
 276     # endif
 277     forAll(edgeFaces, edgeI)
 278     {
 279         //- avoid edges at inter-processor boundaries
 280         if( edgeFaces.sizeOfRow(edgeI) != 2 )
 281             continue;
 282
 283         const label f0 = edgeFaces(edgeI, 0);
 284         const label f1 = edgeFaces(edgeI, 1);
 285
 286         //- both faces have to be part of the same cell
 287         if( faceOwner[f0] != faceOwner[f1] )
 288             continue;
 289
 290         //- check if the feature edge is attachd to exittin faces
 291         if
 292         (
 293             (isBndLayerBase_[f0] && (bFaces[f1].size() == 4)) &&
 294             (isBndLayerBase_[f1] && (bFaces[f0].size() == 4))
 295         )
 296         {
 297             isExitFace_[f0] = true;
 298             isExitFace_[f1] = true;
 299         }
 300     }
 301
 302     # ifdef DEBUGLayer
 303     const label exittingFaceId = mesh_.addFaceSubset("exittingFaces");
 304     forAll(isExitFace_, bfI)
 305         if( isExitFace_[bfI] )
 306             mesh_.addFaceToSubset(exittingFaceId, startBndFaceI+bfI);
 307     # endif
 308
 309     //- classify hair edges as they require different tretment
 310     //- in the smoothing procedure
 311     hairEdgeType_.setSize(hairEdges_.size());
 312
 313     const labelHashSet& corners = mPart.corners();
 314     const labelHashSet& edgePoints = mPart.edgePoints();
 315     const labelHashSet& featureEdges = mPart.featureEdges();
 316
 317     # ifdef USE_OMP
 318     # pragma omp parallel for schedule(dynamic, 100)
 319     # endif
 320     forAll(hairEdges_, hairEdgeI)
 321     {
 322         hairEdgeType_[hairEdgeI] = NONE;
 323
 324         const edge& e = hairEdges_[hairEdgeI];
 325         const label bpI = bp[e.start()];
 326
 327         //- check if it is a boundary edge
 328         forAllRow(bpEdges, bpI, peI)
 329         {
 330             const label beI = bpEdges(bpI, peI);
 331
 332             const edge& be = edges[bpEdges(bpI, peI)];
 333
 334             if( be == e )
 335             {
 336                 hairEdgeType_[hairEdgeI] |= BOUNDARY;
 337
 338                 if( featureEdges.found(beI) )
 339                     hairEdgeType_[hairEdgeI] |= FEATUREEDGE;
 340             }
 341
 342             if( corners.found(bpI) )
 343             {
 344                 hairEdgeType_[hairEdgeI] |= ATCORNER;
 345             }
 346             else if( edgePoints.found(bpI) )
 347             {
 348                 hairEdgeType_[hairEdgeI] |= ATEDGE;
 349             }
 350         }
 351
 352         if( !(hairEdgeType_[hairEdgeI] & BOUNDARY) )
 353             hairEdgeType_[hairEdgeI] |= INSIDE;
 354     }
 355
 356     thinnedHairEdge_.setSize(hairEdges_.size());
 357
 358     //- calculate which other hair edges influence a hair edges
 359     //- and store it in a graph
 360     hairEdgesNearHairEdge_.setSize(hairEdges_.size());
 361
 362     const cellListPMG& cells = mesh_.cells();
 363     const faceList& faces = mesh_.faces();
 364
 365     VRWGraph bpFacesHelper(bpEdges.size());
 366     forAll(faceOwner, bfI)
 367     {
 368         const label cellI = faceOwner[bfI];
 369
 370         const cell& c = cells[cellI];
 371
 372         forAll(c, fI)
 373         {
 374             const face& f = faces[c[fI]];
 375
 376             forAll(f, pI)
 377             {
 378                 const label bpI = bp[f[pI]];
 379                 if( bpI < 0 )
 380                     continue;
 381
 382                 bpFacesHelper.appendIfNotIn(bpI, c[fI]);
 383             }
 384         }
 385     }
 386
 387     forAll(hairEdges_, hairEdgeI)
 388     {
 389         const edge& e = hairEdges_[hairEdgeI];
 390         const label bpI = bp[e.start()];
 391
 392         DynList<label> neiHairEdges;
 393
 394         //- find mesh faces comprising of the current hair edge
 395         forAllRow(bpFacesHelper, bpI, pfI)
 396         {
 397             const face& f = faces[bpFacesHelper(bpI, pfI)];
 398
 399             //- face must be a quad
 400             if( f.size() != 4 )
 401                 continue;
 402
 403             //- check if the current face comprises of the hair edge
 404             label faceEdge(-1);
 405             forAll(f, eI)
 406                 if( f.faceEdge(eI) == e )
 407                 {
 408                     faceEdge = eI;
 409                     break;
 410                 }
 411
 412             if( faceEdge != -1 )
 413             {
 414                 //- check if the opposite edge is also a hair edge
 415                 const label eJ = (faceEdge+2) % 4;
 416
 417                 const edge fe = f.faceEdge(eJ);
 418
 419                 for(label i=0;i<2;++i)
 420                 {
 421                     const label bpJ = bp[fe[i]];
 422
 423                     if( bpJ >= 0 )
 424                     {
 425                         forAllRow(hairEdgesAtBndPoint_, bpJ, pI)
 426                         {
 427                             const label heJ = hairEdgesAtBndPoint_(bpJ, pI);
 428                             if( hairEdges_[heJ] == fe )
 429                                 neiHairEdges.append(heJ);
 430                         }
 431                     }
 432                 }
 433             }
 434         }
 435
 436         hairEdgesNearHairEdge_.setRow(hairEdgeI, neiHairEdges);
 437     }
 438
 439     # ifdef DEBUGLayer
 440     const label hairEdgesId = mesh_.addPointSubset("hairEdgePoints");
 441     const label bndHairEdgeId = mesh_.addPointSubset("bndHairEdgePoints");
 442     const label featureHairEdgeId = mesh_.addPointSubset("featureEdgePoints");
 443     const label cornerHairEdgeId = mesh_.addPointSubset("cornerHairEdgePoints");
 444     const label hairEdgeAtEdgeId = mesh_.addPointSubset("hairEdgeAtEdgePoints");
 445
 446     forAll(hairEdgeType_, heI)
 447     {
 448         const edge& e = hairEdges_[heI];
 449
 450         mesh_.addPointToSubset(hairEdgesId, e.start());
 451         mesh_.addPointToSubset(hairEdgesId, e.end());
 452         if( hairEdgeType_[heI] & FEATUREEDGE)
 453         {
 454             mesh_.addPointToSubset(featureHairEdgeId, e.start());
 455             mesh_.addPointToSubset(featureHairEdgeId, e.end());
 456         }
 457
 458         if( hairEdgeType_[heI] & BOUNDARY)
 459         {
 460             mesh_.addPointToSubset(bndHairEdgeId, e.start());
 461             mesh_.addPointToSubset(bndHairEdgeId, e.end());
 462         }
 463
 464         if( hairEdgeType_[heI] & ATCORNER)
 465         {
 466             mesh_.addPointToSubset(cornerHairEdgeId, e.start());
 467             mesh_.addPointToSubset(cornerHairEdgeId, e.end());
 468         }
 469
 470         if( hairEdgeType_[heI] & ATEDGE)
 471         {
 472             mesh_.addPointToSubset(hairEdgeAtEdgeId, e.start());
 473             mesh_.addPointToSubset(hairEdgeAtEdgeId, e.end());
 474         }
 475     }
 476
 477     mesh_.write();
 478     # endif
 479 }
 480
 481 bool boundaryLayerOptimisation::optimiseLayersAtExittingFaces()
 482 {
 483     bool modified(false);
 484
 485     //- find edge points inside the mesh with more than one hair edge
 486     //- attached to it
 487     labelList nEdgesAtPoint(mesh_.points().size(), 0);
 488
 489     # ifdef USE_OMP
 490     # pragma omp parallel for schedule(dynamic, 50)
 491     # endif
 492     forAll(hairEdges_, heI)
 493     {
 494         # ifdef USE_OMP
 495         # pragma omp atomic
 496         # endif
 497         ++nEdgesAtPoint[hairEdges_[heI].end()];
 498     }
 499
 500     //- find the points with more than one hair edge which was modified
 501     //- in the previous procedure
 502     boolList thinnedPoints(mesh_.points().size(), false);
 503
 504     # ifdef USE_OMP
 505     # pragma omp parallel for schedule(dynamic, 50)
 506     # endif
 507     forAll(thinnedHairEdge_, heI)
 508     {
 509         if
 510         (
 511             thinnedHairEdge_[heI] &&
 512             (nEdgesAtPoint[hairEdges_[heI].end()] > 2)
 513         )
 514         {
 515             modified = true;
 516             thinnedPoints[hairEdges_[heI].end()] = true;
 517         }
 518     }
 519
 520     reduce(modified, maxOp<bool>());
 521
 522     if( !modified )
 523         return false;
 524
 525     Info << "Hair edges at exitting faces shall "
 526          << "be modified due to inner constraints" << endl;
 527
 528     return true;
 529 }
 530
 531 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 532
 533 void boundaryLayerOptimisation::optimiseLayer()
 534 {
 535     //- create surface smoother
 536     meshSurfaceOptimizer surfOpt(meshSurface());
 537
 538     //- lock exitting faces and feature edges
 539     labelLongList lockedFaces;
 540     forAll(isExitFace_, bfI)
 541         if( isExitFace_[bfI] )
 542             lockedFaces.append(bfI);
 543     surfOpt.lockBoundaryFaces(lockedFaces);
 544     surfOpt.lockFeatureEdges();
 545
 546     label nIter(0);
 547     do
 548     {
 549         thinnedHairEdge_ = false;
 550
 551         //- calculate normals at the boundary
 552         optimiseHairNormalsAtTheBoundary();
 553
 554         //- smoothing thickness variation of boundary hairs
 555         optimiseThicknessVariation(BOUNDARY);
 556
 557         if( true )
 558         {
 559             meshSurfaceEngineModifier bMod(meshSurface());
 560             bMod.updateGeometry();
 561
 562             surfOpt.optimizeSurface(2);
 563             bMod.updateGeometry();
 564         }
 565
 566         # ifdef DEBUGLayer
 567         label counter(0);
 568         forAll(thinnedHairEdge_, heI)
 569             if( thinnedHairEdge_[heI] )
 570                 ++counter;
 571         reduce(counter, sumOp<label>());
 572         Info << "Thinned " << counter << " bnd hair edges" << endl;
 573         # endif
 574
 575         //- optimise normals inside the mesh
 576         optimiseHairNormalsInside();
 577
 578         //- optimise thickness variation inside the mesh
 579         optimiseThicknessVariation(INSIDE);
 580
 581         # ifdef DEBUGLayer
 582         label intCounter = 0;
 583         forAll(thinnedHairEdge_, heI)
 584             if( thinnedHairEdge_[heI] )
 585                 ++intCounter;
 586         Info << "Thinned " << (intCounter - counter)
 587              << " inner hair edges" << endl;
 588         # endif
 589     } while( optimiseLayersAtExittingFaces() && (++nIter < maxNumIterations_) );
 590 }
 591
 592 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 593
 594 } // End namespace Foam
 595
 596 // ************************************************************************* //