src/mesh/cfMesh/utilities/smoothers/geometry/meshSurfaceOptimizer/advancedSurfaceSmoothers/surfaceOptimizer/surfaceOptimizer.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 "surfaceOptimizer.H"
  30 #include "matrix2D.H"
  31
  32 //#define DEBUGSmooth
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 const vector surfaceOptimizer::dirVecs[4] =
  42     {
  43         vector(-1.0, -1.0, 0.0),
  44         vector(1.0, -1.0, 0.0),
  45         vector(-1.0, 1.0, 0.0),
  46         vector(1.0, 1.0, 0.0)
  47     };
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 scalar surfaceOptimizer::evaluateStabilisationFactor() const
  52 {
  53     //- find the minimum area
  54     scalar Amin(VGREAT), LsqMax(0.0);
  55     forAll(trias_, triI)
  56     {
  57         const point& p0 = pts_[trias_[triI][0]];
  58         const point& p1 = pts_[trias_[triI][1]];
  59         const point& p2 = pts_[trias_[triI][2]];
  60
  61         const scalar Atri =
  62             0.5 *
  63             (
  64                 (p1.x() - p0.x()) * (p2.y() - p0.y()) -
  65                 (p2.x() - p0.x()) * (p1.y() - p0.y())
  66             );
  67
  68         # ifdef DEBUGSmooth
  69         Info << "Triangle " << triI << " area " << Atri << endl;
  70         # endif
  71
  72         const scalar LSqrTri
  73         (
  74             magSqr(p0 - p1) +
  75             magSqr(p2 - p0)
  76         );
  77
  78         Amin = Foam::min(Amin, Atri);
  79         LsqMax = Foam::max(LsqMax, LSqrTri);
  80     }
  81
  82     # ifdef DEBUGSmooth
  83     Info << "Amin " << Amin << endl;
  84     Info << "LsqMax " << LsqMax << endl;
  85     # endif
  86
  87     //- K is greater than zero in case the stabilisation is needed
  88     scalar K = 0.0;
  89     if( Amin < SMALL * LsqMax )
  90     {
  91         K = SMALL * LsqMax;
  92     }
  93
  94     return K;
  95 }
  96
  97 scalar surfaceOptimizer::evaluateFunc(const scalar& K) const
  98 {
  99     scalar func(0.0);
 100
 101     forAll(trias_, triI)
 102     {
 103         const point& p0 = pts_[trias_[triI][0]];
 104         const point& p1 = pts_[trias_[triI][1]];
 105         const point& p2 = pts_[trias_[triI][2]];
 106
 107         const scalar Atri =
 108             0.5 *
 109             (
 110                 (p1.x() - p0.x()) * (p2.y() - p0.y()) -
 111                 (p2.x() - p0.x()) * (p1.y() - p0.y())
 112             );
 113
 114         const scalar stab = sqrt(sqr(Atri) + K);
 115
 116         # ifdef DEBUGSmooth
 117         Info << "Triangle " << triI << " area " << Atri << endl;
 118         # endif
 119
 120         const scalar LSqrTri
 121         (
 122             magSqr(p0 - p1) +
 123             magSqr(p2 - p0)
 124         );
 125
 126         const scalar Astab = Foam::max(VSMALL, 0.5 * (Atri + stab));
 127         func += LSqrTri / Astab;
 128     }
 129
 130     return func;
 131 }
 132
 133 //- evaluate gradients needed for optimisation
 134 void surfaceOptimizer::evaluateGradients
 135 (
 136     const scalar& K,
 137     vector& gradF,
 138     tensor& gradGradF
 139 ) const
 140 {
 141     gradF = vector::zero;
 142     gradGradF = tensor::zero;
 143
 144     tensor gradGradLt(tensor::zero);
 145     gradGradLt.xx() = 4.0;
 146     gradGradLt.yy() = 4.0;
 147
 148     forAll(trias_, triI)
 149     {
 150         const point& p0 = pts_[trias_[triI][0]];
 151         const point& p1 = pts_[trias_[triI][1]];
 152         const point& p2 = pts_[trias_[triI][2]];
 153
 154         if( magSqr(p1 - p2) < VSMALL ) continue;
 155
 156         const scalar LSqrTri
 157         (
 158             magSqr(p0 - p1) +
 159             magSqr(p2 - p0)
 160         );
 161
 162         const scalar Atri =
 163             0.5 *
 164             (
 165                 (p1.x() - p0.x()) * (p2.y() - p0.y()) -
 166                 (p2.x() - p0.x()) * (p1.y() - p0.y())
 167             );
 168
 169         const scalar stab = sqrt(sqr(Atri) + K);
 170
 171         const scalar Astab = Foam::max(ROOTVSMALL, 0.5 * (Atri + stab));
 172
 173         const vector gradAtri
 174         (
 175             0.5 * (p1.y() - p2.y()),
 176             0.5 * (p2.x() - p1.x()),
 177             0.0
 178         );
 179
 180         const vector gradAstab = 0.5 * (gradAtri + Atri * gradAtri / stab);
 181
 182         const tensor gradGradAstab =
 183             0.5 *
 184             (
 185                 (gradAtri * gradAtri) / stab -
 186                 sqr(Atri) * (gradAtri * gradAtri) / pow(stab, 3)
 187             );
 188
 189         const vector gradLt(4.0 * p0 - 2.0 * p1 - 2.0 * p2);
 190
 191         //- calculate the gradient
 192         const scalar sqrAstab = sqr(Astab);
 193         gradF += gradLt / Astab - (LSqrTri * gradAstab) / sqrAstab;
 194
 195         //- calculate the second gradient
 196         gradGradF +=
 197             gradGradLt / Astab -
 198             twoSymm(gradLt * gradAstab) / sqrAstab -
 199             gradGradAstab * LSqrTri / sqrAstab +
 200             2.0 * LSqrTri * (gradAstab * gradAstab) / (sqrAstab * Astab);
 201     }
 202
 203     //- stabilise diagonal terms
 204     if( mag(gradGradF.xx()) < VSMALL ) gradGradF.xx() = VSMALL;
 205     if( mag(gradGradF.yy()) < VSMALL ) gradGradF.yy() = VSMALL;
 206 }
 207
 208 scalar surfaceOptimizer::optimiseDivideAndConquer(const scalar tol)
 209 {
 210     point& pOpt = pts_[trias_[0][0]];
 211
 212     pOpt = 0.5 * (pMax_ + pMin_);
 213     point currCentre = pOpt;
 214     scalar dx = (pMax_.x() - pMin_.x()) / 2.0;
 215     scalar dy = (pMax_.y() - pMin_.y()) / 2.0;
 216
 217     label iter(0);
 218
 219     //- find the value of the functional in the centre of the bnd box
 220     scalar K = evaluateStabilisationFactor();
 221     scalar funcBefore, funcAfter(evaluateFunc(K));
 222
 223     do
 224     {
 225         funcBefore = funcAfter;
 226
 227         funcAfter = VGREAT;
 228         point minCentre(vector::zero);
 229
 230         for(label i=0;i<4;++i)
 231         {
 232             pOpt.x() = currCentre.x() + 0.5 * dirVecs[i].x() * dx;
 233             pOpt.y() = currCentre.y() + 0.5 * dirVecs[i].y() * dy;
 234
 235             K = evaluateStabilisationFactor();
 236             const scalar func = evaluateFunc(K);
 237
 238             if( func < funcAfter )
 239             {
 240                 minCentre = pOpt;
 241                 funcAfter = func;
 242             }
 243         }
 244
 245         //- set the centre with the minimum value
 246         //- as the centre for future search
 247         currCentre = minCentre;
 248         pOpt = minCentre;
 249
 250         //- halve the search range
 251         dx *= 0.5;
 252         dy *= 0.5;
 253
 254         //- calculate the tolerence
 255         const scalar t = mag(funcAfter - funcBefore) / funcAfter;
 256
 257         # ifdef DEBUGSmooth
 258         Info << "Point position " << pOpt << endl;
 259         Info << "Func before " << funcBefore << endl;
 260         Info << "Func after " << funcAfter << endl;
 261         Info << "Normalised difference " << t << endl;
 262         # endif
 263
 264         if( t < tol )
 265             break;
 266     } while( ++iter < 100 );
 267
 268     return funcAfter;
 269 }
 270
 271 scalar surfaceOptimizer::optimiseSteepestDescent(const scalar tol)
 272 {
 273     point& pOpt = pts_[trias_[0][0]];
 274
 275     //- find the bounding box
 276     const scalar avgEdge = Foam::mag(pMax_ - pMin_);
 277
 278     //- find the minimum value on the 5 x 5 raster
 279     scalar K = evaluateStabilisationFactor();
 280     scalar funcBefore, funcAfter(evaluateFunc(K));
 281
 282     //- start with steepest descent optimisation
 283     vector gradF;
 284     tensor gradGradF;
 285     vector disp;
 286     disp.z() = 0.0;
 287
 288     direction nIterations(0);
 289     do
 290     {
 291         funcBefore = funcAfter;
 292
 293         evaluateGradients(K, gradF, gradGradF);
 294
 295         //- store data into a matrix
 296         matrix2D mat;
 297         mat[0][0] = gradGradF.xx();
 298         mat[0][1] = gradGradF.xy();
 299         mat[1][0] = gradGradF.yx();
 300         mat[1][1] = gradGradF.yy();
 301         FixedList<scalar, 2> source;
 302         source[0] = gradF.x();
 303         source[1] = gradF.y();
 304
 305         //- calculate the determinant
 306         const scalar det = mat.determinant();
 307
 308         if( mag(det) < VSMALL )
 309         {
 310             disp = vector::zero;
 311         }
 312         else
 313         {
 314             disp.x() = mat.solveFirst(source);
 315             disp.y() = mat.solveSecond(source);
 316
 317             if( mag(disp) > 0.2 * avgEdge )
 318             {
 319                 vector dir = disp / mag(disp);
 320
 321                 disp = dir * 0.2 * avgEdge;
 322             }
 323         }
 324
 325         # ifdef DEBUGSmooth
 326         Info << "Second gradient " << gradGradF << endl;
 327         Info << "Gradient " << gradF << endl;
 328         Info << "Displacement " << disp << endl;
 329         Info << "K = " << K << endl;
 330         # endif
 331
 332         pOpt -= disp;
 333
 334         K = evaluateStabilisationFactor();
 335         funcAfter = evaluateFunc(K);
 336
 337         if( mag(funcAfter - funcBefore) / funcBefore < tol )
 338             break;
 339
 340         #ifdef DEBUGSmooth
 341         Info << "New coordinates " << pOpt << endl;
 342         # endif
 343
 344     } while( ++nIterations < 100 );
 345
 346     return funcAfter;
 347 }
 348
 349 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 350
 351 surfaceOptimizer::surfaceOptimizer
 352 (
 353     DynList<point>& pts,
 354     const DynList<triFace>& trias
 355 )
 356 :
 357     pts_(pts),
 358     trias_(trias),
 359     pMin_(),
 360     pMax_()
 361 {
 362     pMin_ = pts_[trias_[0][1]];
 363     pMax_ = pMin_;
 364
 365     forAll(trias_, triI)
 366     {
 367         const triFace& tf = trias_[triI];
 368
 369         for(label i=1;i<3;++i)
 370         {
 371             pMin_ = Foam::min(pMin_, pts_[tf[i]]);
 372             pMax_ = Foam::max(pMax_, pts_[tf[i]]);
 373         }
 374     }
 375 }
 376
 377 surfaceOptimizer::~surfaceOptimizer()
 378 {}
 379
 380 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 381
 382 point surfaceOptimizer::optimizePoint(const scalar tol)
 383 {
 384     const scalar scale = mag(pMax_ - pMin_);
 385     forAll(pts_, i)
 386         pts_[i] /= scale;
 387     pMin_ /= scale;
 388     pMax_ /= scale;
 389
 390     point& pOpt = pts_[trias_[0][0]];
 391
 392     const scalar funcDivide = optimiseDivideAndConquer(tol);
 393     const point newPoint = pOpt;
 394
 395     const scalar funcSteepest = optimiseSteepestDescent(tol);
 396
 397     if( funcSteepest > funcDivide )
 398         pOpt = newPoint;
 399
 400     forAll(pts_, i)
 401         pts_[i] *= scale;
 402     pMin_ *= scale;
 403     pMax_ *= scale;
 404
 405     return pOpt;
 406 }
 407
 408 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 409
 410 } // End namespace Foam
 411
 412 // ************************************************************************* //