src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/boundaryLayerOptimisation/boundaryLayerOptimisationThickness.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 "meshSurfaceEngine.H"
  31 #include "helperFunctions.H"
  32 #include "labelledScalar.H"
  33 #include "polyMeshGenAddressing.H"
  34
  35 //#define DEBUGLayer
  36
  37 # ifdef USE_OMP
  38 #include <omp.h>
  39 # endif
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 void boundaryLayerOptimisation::hairEdgesAtBndFace
  49 (
  50     const label cellI,
  51     const label baseFaceI,
  52     DynList<edge>& hairEdges
  53 ) const
  54 {
  55     const faceListPMG& faces = mesh_.faces();
  56
  57     const cell& c = mesh_.cells()[cellI];
  58
  59         //- check cell topology
  60     DynList<edge, 48> edges;
  61     DynList<DynList<label, 2>, 48> edgeFaces;
  62     DynList<DynList<label, 10>, 24> faceEdges;
  63     faceEdges.setSize(c.size());
  64     label baseFace(-1);
  65     forAll(c, fI)
  66     {
  67         if( c[fI] == baseFaceI )
  68         {
  69             baseFace = fI;
  70         }
  71
  72         const face& f = faces[c[fI]];
  73         faceEdges[fI].setSize(f.size());
  74
  75         forAll(f, eI)
  76         {
  77             const edge e = f.faceEdge(eI);
  78
  79             label pos = edges.containsAtPosition(e);
  80
  81             if( pos < 0 )
  82             {
  83                 pos = edges.size();
  84                 edges.append(e);
  85                 edgeFaces.setSize(pos+1);
  86             }
  87
  88             edgeFaces[pos].append(fI);
  89             faceEdges[fI][eI] = pos;
  90         }
  91     }
  92
  93     const face& bf = faces[c[baseFace]];
  94     hairEdges.setSize(bf.size());
  95
  96     forAll(bf, pI)
  97     {
  98         const label nextEdge = faceEdges[baseFace][pI];
  99         const label prevEdge = faceEdges[baseFace][bf.rcIndex(pI)];
 100
 101         if( edgeFaces[nextEdge].size() != 2 || edgeFaces[prevEdge].size() != 2 )
 102             break;
 103
 104         //- find the face attached to the edge after the current point
 105         label otherNextFace = edgeFaces[nextEdge][0];
 106         if( otherNextFace == baseFace )
 107             otherNextFace = edgeFaces[nextEdge][1];
 108
 109         //- find the face attached to the edge before the current point
 110         label otherPrevFace = edgeFaces[prevEdge][0];
 111         if( otherPrevFace == baseFace )
 112             otherPrevFace = edgeFaces[prevEdge][1];
 113
 114         label commonEdge;
 115         for(commonEdge=0;commonEdge<edges.size();++commonEdge)
 116             if(
 117                 edgeFaces[commonEdge].contains(otherNextFace) &&
 118                 edgeFaces[commonEdge].contains(otherPrevFace)
 119             )
 120                 break;
 121
 122         if( commonEdge == edges.size() )
 123             break;
 124
 125         //- there exists a common edge which shall be used as a hair
 126         if( edges[commonEdge].start() == bf[pI] )
 127         {
 128             hairEdges[pI] = edges[commonEdge];
 129         }
 130         else
 131         {
 132             hairEdges[pI] = edges[commonEdge].reverseEdge();
 133         }
 134     }
 135 }
 136
 137 scalar boundaryLayerOptimisation::calculateThickness
 138 (
 139     const label heI,
 140     const label heJ
 141 ) const
 142 {
 143     const pointFieldPMG& points = mesh_.points();
 144
 145     //- references to hair edges
 146     const edge& he = hairEdges_[heI];
 147     const edge& nhe = hairEdges_[heJ];
 148
 149     //- distance vector between the surface points of hair edges
 150     const point& sp = points[he[0]];
 151     const point& ep = points[nhe[0]];
 152     const vector dv = ep - sp;
 153     const scalar magDv = mag(dv);
 154
 155     //- calculate layer thickness
 156     const scalar currThickness = he.mag(points);
 157     scalar retThickness = currThickness;
 158
 159     const scalar currNeiThickness = nhe.mag(points);
 160     scalar suggestedNeiThickness = currNeiThickness;
 161
 162     //- calculate layer height at the current point
 163     const point npAlpha = help::nearestPointOnTheEdge(sp, ep, points[he[1]]);
 164     const scalar currHeight = mag(npAlpha - points[he[1]]);
 165     scalar retHeight = currHeight;
 166     const scalar cosAlpha = sign((npAlpha - sp) & dv) * mag(npAlpha - sp);
 167     const scalar alpha =
 168         Foam::acos
 169         (
 170             Foam::max
 171             (
 172                 -1.0,
 173                 Foam::min(1.0, cosAlpha / (currThickness + VSMALL))
 174             )
 175         );
 176
 177     //- calculate the height of the layer at the neighbour
 178     //- point
 179     const point npBeta = help::nearestPointOnTheEdge(ep, sp, points[nhe[1]]);
 180     const scalar currNeiHeight = mag(npBeta - points[nhe[1]]);
 181     scalar suggestedNeiHeight = currNeiHeight;
 182     const scalar cosBeta = sign((npBeta - ep) & -dv) * mag(npBeta - ep);
 183     const scalar beta =
 184         Foam::acos
 185         (
 186             Foam::max
 187             (
 188                 -1.0,
 189                 Foam::min(1.0, cosBeta / (currNeiThickness + VSMALL))
 190             )
 191         );
 192
 193     //- check if the current thickness is Ok for the local curvature
 194     if( (alpha + beta) < M_PI )
 195     {
 196         const scalar gamma = M_PI - (alpha + beta);
 197         const scalar sinGamma = Foam::max(SMALL, Foam::sin(gamma));
 198         const scalar sinAlpha = Foam::max(SMALL, Foam::sin(alpha));
 199         const scalar sinBeta = Foam::max(SMALL, Foam::sin(beta));
 200
 201         //- max allowed thickness and layer height due to curvature
 202         retThickness =
 203             Foam::min
 204             (
 205                 retThickness,
 206                 featureSizeFactor_ * magDv * sinBeta / sinGamma
 207             );
 208         retHeight *= (retThickness / (currThickness + VSMALL));
 209
 210         //- max allowed neighbour hair thickness
 211         //- and layer height due to curvature
 212         suggestedNeiThickness =
 213             Foam::min
 214             (
 215                 suggestedNeiThickness,
 216                 featureSizeFactor_ * magDv * sinAlpha / sinGamma
 217             );
 218         suggestedNeiHeight *=
 219             (suggestedNeiThickness / (currNeiThickness + VSMALL));
 220     }
 221
 222     //- check the height variation
 223     const scalar tanVal = (retHeight - suggestedNeiHeight) / (magDv + VSMALL);
 224
 225     if( tanVal > relThicknessTol_ )
 226     {
 227         retHeight = suggestedNeiHeight + relThicknessTol_ * magDv;
 228
 229         retThickness = (retHeight / currHeight) * currThickness;
 230     }
 231
 232     return retThickness;
 233 }
 234
 235 scalar boundaryLayerOptimisation::calculateThicknessOverCell
 236 (
 237     const label heI,
 238     const label cellI,
 239     const label baseFaceI
 240 ) const
 241 {
 242     const pointFieldPMG& points = mesh_.points();
 243     const faceListPMG& faces = mesh_.faces();
 244
 245     const cell& c = mesh_.cells()[cellI];
 246
 247     const face& bf = faces[baseFaceI];
 248
 249     const edge& he = hairEdges_[heI];
 250
 251     const point& sp = points[he[0]];
 252     const point& ep = points[he[1]];
 253
 254     scalar maxThickness = he.mag(points);
 255
 256     //- the base face must not contain the hair edge
 257     //- this is the case at exitting layers
 258     forAll(bf, eI)
 259         if( bf.faceEdge(eI) == he )
 260             return maxThickness;
 261
 262     forAll(c, fI)
 263     {
 264         if( c[fI] == baseFaceI )
 265             continue;
 266
 267         const face& f = faces[c[fI]];
 268
 269         if( help::shareAnEdge(bf, f) && (f.which(he.start()) == -1) )
 270         {
 271             point intersection;
 272
 273             if( !help::lineFaceIntersection(sp, ep, f, points, intersection) )
 274                 continue;
 275
 276             const scalar maxDist = featureSizeFactor_ * mag(intersection - sp);
 277
 278             maxThickness =
 279                 Foam::min(maxThickness, maxDist);
 280         }
 281     }
 282
 283     return maxThickness;
 284 }
 285
 286 void boundaryLayerOptimisation::optimiseThicknessVariation
 287 (
 288     const direction edgeType
 289 )
 290 {
 291     pointFieldPMG& points = mesh_.points();
 292
 293     const meshSurfaceEngine& mse = meshSurface();
 294     const labelList& bp = mse.bp();
 295     const VRWGraph& pFaces = mse.pointFaces();
 296     const faceList::subList& bFaces = mse.boundaryFaces();
 297     const label start = mesh_.nInternalFaces();
 298     const labelList& faceOwner = mse.faceOwners();
 299
 300     vectorField hairDirections(hairEdges_.size());
 301     scalarField hairLength(hairEdges_.size());
 302
 303     # ifdef USE_OMP
 304     # pragma omp parallel for schedule(dynamic, 50)
 305     # endif
 306     forAll(hairEdges_, hairEdgeI)
 307     {
 308         vector n = hairEdges_[hairEdgeI].vec(points);
 309
 310         hairLength[hairEdgeI] = (Foam::mag(n) + VSMALL);
 311         hairDirections[hairEdgeI] = n / hairLength[hairEdgeI];
 312     }
 313
 314     //- reduce thickness of the layer
 315     //- such that the variation of layer thickness
 316     //- It is an iterative process where the layer is thinned in the regions
 317     //- where the tangent is greater than the tolerance value or the curvature
 318     //- permits thicker boundary layers.
 319     boolList activeHairEdge(hairEdges_.size(), true);
 320     bool changed;
 321     label nIter(0);
 322     do
 323     {
 324         changed = false;
 325
 326         boolList modifiedHairEdge(hairEdges_.size(), false);
 327         boolList influencedEdges(hairEdges_.size(), false);
 328
 329         //- check if the hair edge intersects some other face in the cells
 330         //- attached to the hair edge
 331         # ifdef USE_OMP
 332         # pragma omp parallel for schedule(dynamic, 50)
 333         # endif
 334         forAll(hairEdges_, hairEdgeI)
 335         {
 336             if
 337             (
 338                 (hairEdgeType_[hairEdgeI] & edgeType) &&
 339                 activeHairEdge[hairEdgeI]
 340             )
 341             {
 342                 const edge& he = hairEdges_[hairEdgeI];
 343
 344                 const label bpI = bp[he.start()];
 345
 346                 scalar maxThickness = hairLength[hairEdgeI];
 347
 348                 DynList<label, 64> influencers;
 349
 350                 forAllRow(pFaces, bpI, pfI)
 351                 {
 352                     const label bfI = pFaces(bpI, pfI);
 353
 354                     const face& bf = bFaces[bfI];
 355                     if( bf.which(he.end()) >= 0 )
 356                         continue;
 357
 358                     const label baseFaceI = start + bfI;
 359                     const label cOwn = faceOwner[bfI];
 360
 361                     //- check if there exist any self-intersections
 362                     maxThickness =
 363                         Foam::min
 364                         (
 365                             maxThickness,
 366                             calculateThicknessOverCell
 367                             (
 368                                 hairEdgeI,
 369                                 cOwn,
 370                                 baseFaceI
 371                             )
 372                         );
 373
 374                     //- check thickness variation over all hair edges
 375                     DynList<edge> hairEdges;
 376                     hairEdgesAtBndFace(faceOwner[bfI], baseFaceI, hairEdges);
 377
 378                     forAll(bf, pI)
 379                     {
 380                         const edge& nhe = hairEdges[pI];
 381
 382                         const label bpJ = bp[nhe.start()];
 383
 384                         forAllRow(hairEdgesAtBndPoint_, bpJ, peJ)
 385                         {
 386                             const label heJ = hairEdgesAtBndPoint_(bpJ, peJ);
 387
 388                             if( hairEdgeI == heJ )
 389                                 continue;
 390
 391                             if( nhe == hairEdges_[heJ] )
 392                             {
 393                                 influencers.append(heJ);
 394
 395                                 const scalar edgeThickness =
 396                                     calculateThickness(hairEdgeI, heJ);
 397
 398                                 maxThickness =
 399                                     Foam::min(maxThickness, edgeThickness);
 400                             }
 401                         }
 402                     }
 403                 }
 404
 405                 if( hairLength[hairEdgeI] > maxThickness )
 406                 {
 407                     //- make the hair edge shorter
 408                     hairLength[hairEdgeI] = maxThickness;
 409                     modifiedHairEdge[hairEdgeI] = true;
 410                     changed = true;
 411
 412                     thinnedHairEdge_[hairEdgeI] = true;
 413
 414                     forAll(influencers, i)
 415                         influencedEdges[influencers[i]] = true;
 416                 }
 417             }
 418         }
 419 /*
 420         # ifdef USE_OMP
 421         # pragma omp parallel for schedule(dynamic, 50)
 422         # endif
 423         forAll(hairEdgesNearHairEdge_, hairEdgeI)
 424         {
 425             const scalar magN = hairLength[hairEdgeI];
 426
 427             if( magN < VSMALL )
 428                 FatalErrorIn
 429                 (
 430                     "void boundaryLayerOptimisation::optimiseThicknessVariation"
 431                     "(const direction, const label, const scalar, const scalar)"
 432                 ) << "Zero layer thickness at hair edge " << hairEdgeI
 433                   << ". Exitting..." << exit(FatalError);
 434
 435             if( hairEdgeType_[hairEdgeI] & edgeType )
 436             {
 437                 forAllRow(hairEdgesNearHairEdge_, hairEdgeI, nheI)
 438                 {
 439                     const label hairEdgeJ =
 440                         hairEdgesNearHairEdge_(hairEdgeI, nheI);
 441
 442                     if( !activeHairEdge[hairEdgeJ] )
 443                         continue;
 444
 445                     const scalar maxThickness =
 446                         calculateThickness
 447                         (
 448                             hairEdgeI,
 449                             hairEdgeJ
 450                         );
 451
 452                     if( hairLength[hairEdgeI] > maxThickness )
 453                     {
 454                         //- make the hair edge shorter
 455                         hairLength[hairEdgeI] = maxThickness;
 456
 457                         changed = true;
 458                         thinnedHairEdge_[hairEdgeI] = true;
 459                         modifiedHairEdge[hairEdgeI] = true;
 460                     }
 461                 }
 462             }
 463         }
 464         */
 465
 466         if( Pstream::parRun() )
 467         {
 468             //- collect data at inter-processor boundaries
 469             const Map<label>& globalToLocal =
 470                 mse.globalToLocalBndPointAddressing();
 471
 472             const edgeList& edges = mesh_.addressingData().edges();
 473             const VRWGraph& pointEdges = mesh_.addressingData().pointEdges();
 474             const VRWGraph& edgesAtProcs =
 475                 mesh_.addressingData().edgeAtProcs();
 476             const labelLongList& globalEdgeLabel =
 477                 mesh_.addressingData().globalEdgeLabel();
 478             const Map<label>& globalToLocalEdge =
 479                 mesh_.addressingData().globalToLocalEdgeAddressing();
 480             const DynList<label>& eNeiProcs =
 481                 mesh_.addressingData().edgeNeiProcs();
 482
 483             std::map<label, LongList<labelledScalar> > exchangeData;
 484             forAll(eNeiProcs, i)
 485                 exchangeData[eNeiProcs[i]].clear();
 486
 487             forAllConstIter(Map<label>, globalToLocal, it)
 488             {
 489                 const label bpI = it();
 490
 491                 forAllRow(hairEdgesAtBndPoint_, bpI, i)
 492                 {
 493                     const label hairEdgeI = hairEdgesAtBndPoint_(bpI, i);
 494
 495                     if( !(hairEdgeType_[hairEdgeI] & edgeType) )
 496                         continue;
 497
 498                     const edge& he = hairEdges_[hairEdgeI];
 499
 500                     forAllRow(pointEdges, he.start(), peI)
 501                     {
 502                         const label edgeI = pointEdges(he.start(), peI);
 503
 504                         if( he == edges[edgeI] )
 505                         {
 506                             labelledScalar lScalar
 507                             (
 508                                 globalEdgeLabel[edgeI],
 509                                 hairLength[hairEdgeI]
 510                             );
 511
 512                             forAllRow(edgesAtProcs, edgeI, j)
 513                             {
 514                                 const label neiProc = edgesAtProcs(edgeI, j);
 515
 516                                 if( neiProc == Pstream::myProcNo() )
 517                                     continue;
 518
 519                                 LongList<labelledScalar>& dts =
 520                                     exchangeData[neiProc];
 521
 522                                 dts.append(lScalar);
 523                             }
 524                         }
 525                     }
 526                 }
 527             }
 528
 529             LongList<labelledScalar> receivedData;
 530             help::exchangeMap(exchangeData, receivedData);
 531
 532             forAll(receivedData, i)
 533             {
 534                 const labelledScalar& lScalar = receivedData[i];
 535                 const label edgeI = globalToLocalEdge[lScalar.scalarLabel()];
 536                 const edge& e = edges[edgeI];
 537
 538                 bool found(false);
 539                 for(label pI=0;pI<2;++pI)
 540                 {
 541                     const label bpI = bp[e[pI]];
 542
 543                     if( bpI < 0 )
 544                         continue;
 545
 546                     forAllRow(hairEdgesAtBndPoint_, bpI, i)
 547                     {
 548                         const label hairEdgeI = hairEdgesAtBndPoint_(bpI, i);
 549
 550                         if( hairEdges_[hairEdgeI] == e )
 551                         {
 552                             if( lScalar.value() < hairLength[hairEdgeI] )
 553                             {
 554                                 hairLength[hairEdgeI] = lScalar.value();
 555                                 changed = true;
 556                                 thinnedHairEdge_[hairEdgeI] = true;
 557                                 modifiedHairEdge[hairEdgeI] = true;
 558                             }
 559
 560                             found = true;
 561                             break;
 562                         }
 563                     }
 564
 565                     if( found )
 566                         break;
 567                 }
 568             }
 569         }
 570
 571         //- reduce the information over all processors
 572         reduce(changed, maxOp<bool>());
 573
 574         if( !changed )
 575             break;
 576
 577         //- move boundary vertices to the new positions
 578         # ifdef USE_OMP
 579         # pragma omp parallel for schedule(dynamic, 100)
 580         # endif
 581         forAll(hairEdges_, hairEdgeI)
 582         {
 583             if( hairEdgeType_[hairEdgeI] & edgeType )
 584             {
 585                 const edge& he = hairEdges_[hairEdgeI];
 586                 const vector& hv = hairDirections[hairEdgeI];
 587                 const point& s = points[he.start()];
 588
 589                 points[he.end()] = s + hairLength[hairEdgeI] * hv;
 590             }
 591         }
 592
 593         //- mark edges which may be changed
 594         activeHairEdge.transfer(influencedEdges);
 595     } while( changed && (++nIter < 1000) );
 596 }
 597
 598
 599 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 600
 601 } // End namespace Foam
 602
 603 // ************************************************************************* //