src/mesh/cfMesh/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerI.H

   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 "meshSurfaceOptimizer.H"
  30 #include "plane.H"
  31 #include "Map.H"
  32 #include "surfaceOptimizer.H"
  33 #include "helperFunctions.H"
  34 #include "partTriMeshSimplex.H"
  35
  36 //#define DEBUGSmooth
  37
  38 # ifdef DEBUGSmooth
  39 #include "triSurf.H"
  40 #include "triSurfModifier.H"
  41 #include "helperFunctions.H"
  42 # endif
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 inline const partTriMesh& meshSurfaceOptimizer::triMesh() const
  52 {
  53     if( !triMeshPtr_ )
  54         calculateTrianglesAndAddressing();
  55
  56     return *triMeshPtr_;
  57 }
  58
  59 inline void meshSurfaceOptimizer::updateTriMesh(const labelLongList& selPoints)
  60 {
  61     if( !triMeshPtr_ )
  62         FatalErrorIn
  63         (
  64             "inline void meshSurfaceOptimizer::updateTriMesh"
  65             "(const labelLongList&)"
  66         ) << "triMeshPtr_ is not allocated " << abort(FatalError);
  67
  68     triMeshPtr_->updateVertices(selPoints);
  69 }
  70
  71 inline void meshSurfaceOptimizer::updateTriMesh()
  72 {
  73     if( !triMeshPtr_ )
  74         FatalErrorIn
  75         (
  76             "inline void meshSurfaceOptimizer::updateTriMesh()"
  77         ) << "triMeshPtr_ is not allocated " << abort(FatalError);
  78
  79     triMeshPtr_->updateVertices();
  80 }
  81
  82 inline bool meshSurfaceOptimizer::transformIntoPlane
  83 (
  84     const label bpI,
  85     const plane& pl,
  86     vector& vecX,
  87     vector& vecY,
  88     DynList<point>& pts,
  89     DynList<triFace>& trias
  90 ) const
  91 {
  92     # ifdef DEBUGSmooth
  93     Pout << "Transforming boundary node " << bpI << endl;
  94     # endif
  95
  96     const partTriMesh& triMesh = this->triMesh();
  97     const label triPointI = triMesh.meshSurfacePointLabelInTriMesh()[bpI];
  98     if( triPointI < 0 )
  99         return false;
 100
 101     partTriMeshSimplex simplex(triMesh, triPointI);
 102
 103     const DynList<point, 32>& sPts = simplex.pts();
 104
 105     //- create vecX and vecY
 106     const point& p = simplex.centrePoint();
 107     pts.setSize(0);
 108     bool found(false);
 109     forAll(sPts, pI)
 110     {
 111         const point sp = pl.nearestPoint(sPts[pI]);
 112         const scalar d = mag(sp - p);
 113         if( d > VSMALL )
 114         {
 115             vecX = sp - p;
 116             vecX /= d;
 117             vecY = pl.normal() ^ vecX;
 118             vecY /= mag(vecY);
 119             found = true;
 120             break;
 121         }
 122     }
 123
 124     if( !found )
 125         return false;
 126
 127     trias = simplex.triangles();
 128
 129     # ifdef DEBUGSmooth
 130     Pout << "VecX " << vecX << endl;
 131     Pout << "vecY " << vecY << endl;
 132     # endif
 133
 134     //- transform the vertices
 135     pts.setSize(sPts.size());
 136
 137     forAll(sPts, pI)
 138     {
 139         const point pt
 140         (
 141             ((sPts[pI] - p) & vecX),
 142             ((sPts[pI] - p) & vecY),
 143             0.0
 144         );
 145
 146         pts[pI] = pt;
 147     }
 148
 149     # ifdef DEBUGSmooth
 150     Info << "Original triangles " << endl;
 151     forAll(simplex.triangles(), triI)
 152         Info << "Tri " << triI << " is " << simplex.triangles()[triI] << endl;
 153     Info << "Transformed triangles are " << trias << endl;
 154     Info << "Transformed vertices " << pts << endl;
 155
 156     triSurf surf;
 157     triSurfModifier sMod(surf);
 158     pointField& sPoints = sMod.pointsAccess();
 159     sPoints.setSize(pts.size());
 160     forAll(sPoints, i)
 161         sPoints[i] = pts[i];
 162     LongList<labelledTri>& sTrias = sMod.facetsAccess();
 163     sTrias.setSize(trias.size());
 164     forAll(sTrias, i)
 165     {
 166         labelledTri& tf = sTrias[i];
 167         tf[0] = trias[i][0];
 168         tf[1] = trias[i][1];
 169         tf[2] = trias[i][2];
 170
 171         tf.region() = 0;
 172     }
 173     sMod.patchesAccess().setSize(1);
 174     sMod.patchesAccess()[0].name() = "bnd";
 175     sMod.patchesAccess()[0].geometricType() = "patch";
 176
 177     fileName sName("bndSimplex_");
 178     sName += help::scalarToText(bpI);
 179     sName += ".stl";
 180     surf.writeSurface(sName);
 181     # endif
 182
 183     return found;
 184 }
 185
 186 inline point meshSurfaceOptimizer::newPositionLaplacian
 187 (
 188     const label bpI,
 189     const bool transformIntoPlane
 190 ) const
 191 {
 192     const VRWGraph& pPoints = surfaceEngine_.pointPoints();
 193     const pointFieldPMG& points = surfaceEngine_.points();
 194     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 195
 196     if( vertexType_[bpI] & LOCKED )
 197         return points[bPoints[bpI]];
 198
 199     vector newP(vector::zero);
 200     if( transformIntoPlane )
 201     {
 202         const vector& pNormal = surfaceEngine_.pointNormals()[bpI];
 203
 204         if( magSqr(pNormal) < VSMALL )
 205             return points[bPoints[bpI]];
 206
 207         plane pl(points[bPoints[bpI]], pNormal);
 208
 209         DynList<point> projectedPoints;
 210         projectedPoints.setSize(pPoints.sizeOfRow(bpI));
 211         forAllRow(pPoints, bpI, pI)
 212             projectedPoints[pI] =
 213                 pl.nearestPoint(points[bPoints[pPoints(bpI, pI)]]);
 214
 215         forAll(projectedPoints, pI)
 216             newP += projectedPoints[pI];
 217
 218         newP /= projectedPoints.size();
 219     }
 220     else
 221     {
 222         forAllRow(pPoints, bpI, pI)
 223             newP += points[bPoints[pPoints(bpI, pI)]];
 224
 225         newP /= pPoints.sizeOfRow(bpI);
 226     }
 227
 228     return newP;
 229 }
 230
 231 inline point meshSurfaceOptimizer::newPositionLaplacianFC
 232 (
 233     const label bpI,
 234     const bool transformIntoPlane
 235 ) const
 236 {
 237     const VRWGraph& pointFaces = surfaceEngine_.pointFaces();
 238     const pointFieldPMG& points = surfaceEngine_.points();
 239     const vectorField& faceCentres = surfaceEngine_.faceCentres();
 240     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 241
 242     if( vertexType_[bpI] & LOCKED )
 243         return points[bPoints[bpI]];
 244
 245     vector newP(vector::zero);
 246
 247     if( transformIntoPlane )
 248     {
 249         const vector& pNormal = surfaceEngine_.pointNormals()[bpI];
 250
 251         if( magSqr(pNormal) < VSMALL )
 252             return points[bPoints[bpI]];
 253
 254         plane pl(points[bPoints[bpI]], pNormal);
 255
 256         DynList<point> projectedPoints;
 257         projectedPoints.setSize(pointFaces.sizeOfRow(bpI));
 258         forAllRow(pointFaces, bpI, pfI)
 259             projectedPoints[pfI] =
 260                 pl.nearestPoint(faceCentres[pointFaces(bpI, pfI)]);
 261
 262         forAll(projectedPoints, pI)
 263             newP += projectedPoints[pI];
 264
 265         newP /= projectedPoints.size();
 266     }
 267     else
 268     {
 269         forAllRow(pointFaces, bpI, pfI)
 270             newP += faceCentres[pointFaces(bpI, pfI)];
 271
 272         newP /= pointFaces.sizeOfRow(bpI);
 273     }
 274
 275     return newP;
 276 }
 277
 278 inline point meshSurfaceOptimizer::newPositionLaplacianWFC
 279 (
 280     const label bpI,
 281     const bool transformIntoPlane
 282 ) const
 283 {
 284     const VRWGraph& pointFaces = surfaceEngine_.pointFaces();
 285     const pointFieldPMG& points = surfaceEngine_.points();
 286     const vectorField& faceAreas = surfaceEngine_.faceNormals();
 287     const vectorField& faceCentres = surfaceEngine_.faceCentres();
 288     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 289
 290     if( vertexType_[bpI] & LOCKED )
 291         return points[bPoints[bpI]];
 292
 293     vector newP(vector::zero);
 294
 295     if( transformIntoPlane )
 296     {
 297         const vector& pNormal = surfaceEngine_.pointNormals()[bpI];
 298
 299         if( magSqr(pNormal) < VSMALL )
 300             return points[bPoints[bpI]];
 301
 302         plane pl(points[bPoints[bpI]], pNormal);
 303
 304         DynList<point> projectedPoints;
 305         projectedPoints.setSize(pointFaces.sizeOfRow(bpI));
 306         forAllRow(pointFaces, bpI, pfI)
 307             projectedPoints[pfI] =
 308                 pl.nearestPoint(faceCentres[pointFaces(bpI, pfI)]);
 309
 310         scalar sumW(0.0);
 311         forAll(projectedPoints, pI)
 312         {
 313             const label bfI = pointFaces(bpI, pI);
 314             const scalar w = (faceAreas[bfI] & faceAreas[bfI]);
 315             newP += w * projectedPoints[pI];
 316             sumW += w;
 317         }
 318
 319         newP /= Foam::max(sumW, VSMALL);
 320     }
 321     else
 322     {
 323         scalar sumW(0.0);
 324         forAllRow(pointFaces, bpI, pfI)
 325         {
 326             const label bfI = pointFaces(bpI, pfI);
 327             const scalar w = (faceAreas[bfI] & faceAreas[bfI]);
 328             newP += w * faceCentres[pointFaces(bpI, pfI)];
 329             sumW += w;
 330         }
 331
 332         newP /= Foam::max(sumW, VSMALL);
 333     }
 334
 335     return newP;
 336 }
 337
 338 inline point meshSurfaceOptimizer::newPositionSurfaceOptimizer
 339 (
 340     const label bpI,
 341     const scalar tol
 342 ) const
 343 {
 344     const pointFieldPMG& points = surfaceEngine_.points();
 345     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 346
 347     if( vertexType_[bpI] & LOCKED )
 348         return points[bPoints[bpI]];
 349
 350     # ifdef DEBUGSmooth
 351     Pout << "Smoothing boundary node " << bpI << endl;
 352     Pout << "Node label in the mesh is " << bPoints[bpI] << endl;
 353     Pout << "Point coordinates " << points[bPoints[bpI]] << endl;
 354     # endif
 355
 356     //- project vertices onto the plane
 357     const vector& pNormal = surfaceEngine_.pointNormals()[bpI];
 358     if( magSqr(pNormal) < VSMALL )
 359         return points[bPoints[bpI]];
 360
 361     const plane pl(points[bPoints[bpI]], pNormal);
 362
 363     DynList<point> pts;
 364     DynList<triFace> trias;
 365     vector vecX, vecY;
 366     bool success = this->transformIntoPlane(bpI, pl, vecX, vecY, pts, trias);
 367     if( !success )
 368     {
 369         //Warning << "Cannot transform to plane" << endl;
 370         return points[bPoints[bpI]];
 371     }
 372
 373     surfaceOptimizer so(pts, trias);
 374     const point newPoint = so.optimizePoint(tol);
 375
 376     const point newP
 377     (
 378         points[bPoints[bpI]] +
 379         vecX * newPoint.x() +
 380         vecY * newPoint.y()
 381     );
 382
 383     if( help::isnan(newP) || help::isinf(newP) )
 384     {
 385         WarningIn
 386         (
 387             "inline point meshSurfaceOptimizer::newPositionSurfaceOptimizer"
 388             "(const label, const scalar) const"
 389         ) << "Cannot move point " << bpI << endl;
 390
 391         return points[bPoints[bpI]];
 392     }
 393
 394     return newP;
 395 }
 396
 397 inline point meshSurfaceOptimizer::newEdgePositionLaplacian
 398 (
 399     const label bpI
 400 ) const
 401 {
 402     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 403     const edgeList& edges = surfaceEngine_.edges();
 404     const VRWGraph& bpEdges = surfaceEngine_.boundaryPointEdges();
 405     const pointFieldPMG& points = surfaceEngine_.points();
 406
 407     if( vertexType_[bpI] & LOCKED )
 408         return points[bPoints[bpI]];
 409
 410     const labelHashSet& featureEdges = partitionerPtr_->featureEdges();
 411
 412     DynList<label> edgePoints;
 413
 414     forAllRow(bpEdges, bpI, i)
 415     {
 416         const label beI = bpEdges(bpI, i);
 417
 418         if( featureEdges.found(beI) )
 419         {
 420             edgePoints.append(edges[beI].otherVertex(bPoints[bpI]));
 421         }
 422     }
 423
 424     if( edgePoints.size() != 2 )
 425         return points[bPoints[bpI]];
 426
 427     # ifdef DEBUGSearch
 428     Info << "Edge points " << edgePoints << endl;
 429     # endif
 430
 431     vector pos(vector::zero);
 432     forAll(edgePoints, epI)
 433         pos += points[edgePoints[epI]];
 434
 435     pos /= edgePoints.size();
 436
 437     return pos;
 438 }
 439
 440 template<class labelListType>
 441 inline void meshSurfaceOptimizer::lockBoundaryFaces(const labelListType& l)
 442 {
 443     lockedSurfaceFaces_ = l;
 444
 445     const faceList::subList& bFaces = surfaceEngine_.boundaryFaces();
 446     const labelList& bp = surfaceEngine_.bp();
 447
 448     # ifdef USE_OMP
 449     # pragma omp parallel for schedule(dynamic, 20)
 450     # endif
 451     forAll(lockedSurfaceFaces_, lfI)
 452     {
 453         const face& bf = bFaces[lockedSurfaceFaces_[lfI]];
 454
 455         forAll(bf, pI)
 456             vertexType_[bp[bf[pI]]] |= LOCKED;
 457     }
 458
 459     if( Pstream::parRun() )
 460     {
 461         const Map<label>& globalToLocal =
 462             surfaceEngine_.globalToLocalBndPointAddressing();
 463         const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 464         const DynList<label>& bpNeiProcs = surfaceEngine_.bpNeiProcs();
 465
 466         std::map<label, labelLongList> exchangeData;
 467         forAll(bpNeiProcs, i)
 468             exchangeData[bpNeiProcs[i]].clear();
 469
 470         //- prepare data which will be sent to other processors
 471         forAllConstIter(Map<label>, globalToLocal, it)
 472         {
 473             const label bpI = it();
 474
 475             if( vertexType_[bpI] & LOCKED )
 476             {
 477                 forAllRow(bpAtProcs, bpI, i)
 478                 {
 479                     const label neiProc = bpAtProcs(bpI, i);
 480
 481                     if( neiProc == Pstream::myProcNo() )
 482                         continue;
 483
 484                     exchangeData[neiProc].append(it.key());
 485                 }
 486             }
 487         }
 488
 489         labelLongList receivedData;
 490         help::exchangeMap(exchangeData, receivedData);
 491
 492         forAll(receivedData, i)
 493         {
 494             const label bpI = globalToLocal[receivedData[i]];
 495
 496             vertexType_[bpI] |= LOCKED;
 497         }
 498     }
 499 }
 500
 501 template<class labelListType>
 502 inline void meshSurfaceOptimizer::lockBoundaryPoints(const labelListType& l)
 503 {
 504     # ifdef USE_OMP
 505     # pragma omp parallel for schedule(dynamic, 50)
 506     # endif
 507     forAll(l, i)
 508     {
 509         const label bpI = l[i];
 510
 511         vertexType_[bpI] |= LOCKED;
 512     }
 513
 514     if( Pstream::parRun() )
 515     {
 516         const Map<label>& globalToLocal =
 517             surfaceEngine_.globalToLocalBndPointAddressing();
 518         const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 519         const DynList<label>& bpNeiProcs = surfaceEngine_.bpNeiProcs();
 520
 521         std::map<label, labelLongList> exchangeData;
 522         forAll(bpNeiProcs, i)
 523             exchangeData[bpNeiProcs[i]].clear();
 524
 525         //- prepare data which will be sent to other processors
 526         forAllConstIter(Map<label>, globalToLocal, it)
 527         {
 528             const label bpI = it();
 529
 530             if( vertexType_[bpI] & LOCKED )
 531             {
 532                 forAllRow(bpAtProcs, bpI, i)
 533                 {
 534                     const label neiProc = bpAtProcs(bpI, i);
 535
 536                     if( neiProc == Pstream::myProcNo() )
 537                         continue;
 538
 539                     exchangeData[neiProc].append(it.key());
 540                 }
 541             }
 542         }
 543
 544         labelLongList receivedData;
 545         help::exchangeMap(exchangeData, receivedData);
 546
 547         forAll(receivedData, i)
 548         {
 549             const label bpI = globalToLocal[receivedData[i]];
 550
 551             vertexType_[bpI] |= LOCKED;
 552         }
 553     }
 554 }
 555
 556 //- lock edge points
 557 inline void meshSurfaceOptimizer::lockFeatureEdges()
 558 {
 559     # ifdef USE_OMP
 560     # pragma omp parallel for schedule(dynamic, 50)
 561     # endif
 562     forAll(vertexType_, bpI)
 563         if( vertexType_[bpI] & (EDGE | CORNER) )
 564             vertexType_[bpI] |= LOCKED;
 565 }
 566
 567 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 568
 569 } // End namespace Foam
 570
 571 // ************************************************************************* //