src/mesh/cfMesh/utilities/triSurfaceTools/triSurfaceCurvatureEstimator/triSurfaceCurvatureEstimatorCalculate.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 "triSurfaceCurvatureEstimator.H"
  29 #include "matrix3D.H"
  30 #include "quadricFitting.H"
  31 #include "helperFunctions.H"
  32 #include "HashSet.H"
  33 #include "boolList.H"
  34 #include "Map.H"
  35 #include "DynList.H"
  36
  37 #include "OFstream.H"
  38
  39 #include <map>
  40
  41 # ifdef USE_OMP
  42 #include <omp.h>
  43 # endif
  44
  45 //#define DEBUGCurvatureEstimator
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 namespace Foam
  50 {
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53
  54 void writeSurfaceToVTK
  55 (
  56     OFstream& file,
  57     const triSurf& surf
  58 )
  59 {
  60     //- write the header
  61     file << "# vtk DataFile Version 3.0\n";
  62     file << "vtk output\n";
  63     file << "ASCII\n";
  64     file << "DATASET POLYDATA\n";
  65
  66     //- write points
  67     const pointField& points = surf.points();
  68     file << "POINTS " << points.size() << " float\n";
  69     forAll(points, pI)
  70     {
  71         const point& p = points[pI];
  72
  73         file << p.x() << ' ' << p.y() << ' ' << p.z() << ' ';
  74
  75         if( pI % 5 == 0 )
  76             file << "\n";
  77     }
  78
  79     //- write triangles
  80     file << "\n";
  81     file << "\nPOLYGONS " << surf.size() << " " << 4*surf.size() << endl;
  82     forAll(surf, triI)
  83     {
  84         const labelledTri& tri = surf[triI];
  85         file << 3 << " " << tri[0] << " " << tri[1] << " " << tri[2] << endl;
  86     }
  87 }
  88
  89 void writeSurfaceToVTK
  90 (
  91     OFstream& file,
  92     const triSurf& surf,
  93     const DynList<label>& triangles,
  94     const labelHashSet& labels
  95 )
  96 {
  97     //- write the header
  98     file << "# vtk DataFile Version 3.0\n";
  99     file << "vtk output\n";
 100     file << "ASCII\n";
 101     file << "DATASET POLYDATA\n";
 102
 103     //- write points
 104     std::map<label, label> newPointLabel;
 105     forAll(triangles, tI)
 106     {
 107         const labelledTri& tri = surf[triangles[tI]];
 108
 109         for(label pI=0;pI<3;++pI)
 110         {
 111             newPointLabel[tri[pI]];
 112         }
 113     }
 114
 115     forAllConstIter(labelHashSet, labels, it)
 116         newPointLabel[it.key()];
 117
 118     const pointField& points = surf.points();
 119     file << "POINTS " << label(newPointLabel.size()) << " float\n";
 120     label counter(0);
 121     for
 122     (
 123         std::map<label, label>::iterator it=newPointLabel.begin();
 124         it!=newPointLabel.end();
 125         ++it
 126     )
 127     {
 128         it->second = counter++;
 129
 130         const point& p = points[it->first];
 131
 132         file << p.x() << ' ' << p.y() << ' ' << p.z() << ' ';
 133
 134         file << nl;
 135     }
 136
 137     //- write triangles
 138     file << "\n";
 139     file << "\nPOLYGONS " << triangles.size()
 140          << " " << 4*triangles.size() << endl;
 141     forAll(triangles, tI)
 142     {
 143         const labelledTri& tri = surf[triangles[tI]];
 144
 145         file << 3
 146              << " " << newPointLabel[tri[0]]
 147              << " " << newPointLabel[tri[1]]
 148              << " " << newPointLabel[tri[2]] << endl;
 149     }
 150 }
 151
 152 void writeSurfaceToVTK
 153 (
 154     const word& name,
 155     const triSurf& surf,
 156     const List<DynList<scalar, 1> >& data
 157 )
 158 {
 159     OFstream file(name+".vtk");
 160
 161     writeSurfaceToVTK(file, surf);
 162
 163     //- write curvature fields
 164     const pointField& points = surf.points();
 165     file << "\n";
 166     file << "\nPOINT_DATA " << points.size() << "\n";
 167
 168     file << "SCALARS Values double\n";
 169     file << "LOOKUP_TABLE default\n";
 170     forAll(points, pI)
 171     {
 172         file << data[pI][0] << " ";
 173
 174         if( pI && pI % 5 == 0 )
 175             file << endl;
 176     }
 177 }
 178
 179 void writeSurfaceToVTK
 180 (
 181     const word& name,
 182     const triSurf& surf,
 183     const List<DynList<vector, 1> >& data
 184 )
 185 {
 186     OFstream file(name+".vtk");
 187
 188     writeSurfaceToVTK(file, surf);
 189
 190     //- write curvature fields
 191     const pointField& points = surf.points();
 192     file << "\n";
 193     file << "\nPOINT_DATA " << points.size() << "\n";
 194
 195     file << "VECTORS Values double\n";
 196     file << "LOOKUP_TABLE default\n";
 197     forAll(points, pI)
 198     {
 199         const vector& v = data[pI][0];
 200
 201         file << v[0] << " " << v[1] << " " << v[2] << " ";
 202
 203         if( pI && pI % 5 == 0 )
 204             file << endl;
 205     }
 206 }
 207
 208 void triSurfaceCurvatureEstimator::calculateEdgeCurvature()
 209 {
 210     const pointField& points = surface_.points();
 211     const edgeLongList& edges = surface_.edges();
 212     const VRWGraph& pointEdges = surface_.pointEdges();
 213     const VRWGraph& edgeFaces = surface_.edgeFacets();
 214
 215     edgePointCurvature_.setSize(points.size());
 216     boolList featureEdge(edges.size());
 217
 218     # ifdef USE_OMP
 219     # pragma omp parallel
 220     # endif
 221     {
 222         # ifdef USE_OMP
 223         # pragma omp for schedule(static, 1)
 224         # endif
 225         forAll(edgePointCurvature_, i)
 226             edgePointCurvature_[i] = 0.0;
 227
 228         //- mark feature edges
 229         # ifdef USE_OMP
 230         # pragma omp for schedule(static, 1)
 231         # endif
 232         forAll(edgeFaces, eI)
 233         {
 234             if( edgeFaces.sizeOfRow(eI) != 2 )
 235             {
 236                 featureEdge[eI] = false;
 237                 continue;
 238             }
 239
 240             if(
 241                  surface_[edgeFaces(eI, 0)].region() ==
 242                  surface_[edgeFaces(eI, 1)].region()
 243             )
 244             {
 245                 featureEdge[eI] = false;
 246             }
 247             else
 248             {
 249                 featureEdge[eI] = true;
 250             }
 251         }
 252
 253         # ifdef USE_OMP
 254         # pragma omp barrier
 255         # endif
 256
 257         //- loop through the points and calculate the curvature for points
 258         //- attached to two feature edges
 259         # ifdef USE_OMP
 260         # pragma omp for schedule(dynamic, 20)
 261         # endif
 262         forAll(pointEdges, pI)
 263         {
 264             DynList<label> features;
 265             forAllRow(pointEdges, pI, peI)
 266             {
 267                 const label edgeI = pointEdges(pI, peI);
 268                 if( featureEdge[edgeI] )
 269                     features.append(edgeI);
 270             }
 271
 272             if( features.size() == 2 )
 273             {
 274                 //- calculate the curvature and store it in the map
 275                 vector e1 = edges[features[0]].vec(points);
 276                 const scalar d1 = mag(e1) + VSMALL;
 277                 e1 /= d1;
 278                 vector e2 = edges[features[1]].vec(points);
 279                 const scalar d2 = mag(e2) + VSMALL;
 280                 e2 /= d2;
 281
 282                 scalar cs = e1 & e2;
 283                 cs = Foam::min(1.0, cs);
 284                 cs = Foam::max(-1.0, cs);
 285
 286                 const scalar curv = Foam::acos(cs) / (0.5 * (d1+d2+VSMALL));
 287                 edgePointCurvature_[pI] = Foam::mag(curv);
 288             }
 289         }
 290     }
 291 }
 292
 293 void triSurfaceCurvatureEstimator::calculateSurfaceCurvatures()
 294 {
 295     const VRWGraph& pointTriangles = surface_.pointFacets();
 296
 297     const pointField& points = surface_.points();
 298
 299     patchPositions_.setSize(surface_.size());
 300     gaussianCurvature_.setSize(points.size());
 301     meanCurvature_.setSize(points.size());
 302     maxCurvature_.setSize(points.size());
 303     minCurvature_.setSize(points.size());
 304     maxCurvatureVector_.setSize(points.size());
 305     minCurvatureVector_.setSize(points.size());
 306
 307     List<DynList<label, 4> > pointPatches(points.size());
 308
 309     # ifdef USE_OMP
 310     # pragma omp parallel for schedule(dynamic, 40)
 311     # endif
 312     forAll(pointTriangles, pointI)
 313     {
 314         std::map<label, DynList<label> > regionTriangles;
 315         Map<labelHashSet> otherLabels;
 316         Map<vector> normals;
 317
 318         forAllRow(pointTriangles, pointI, ptI)
 319         {
 320             const label triI = pointTriangles(pointI, ptI);
 321             const label regionI = surface_[triI].region();
 322
 323             if( !normals.found(regionI) )
 324                 normals.insert(regionI, vector::zero);
 325             if( regionTriangles.find(regionI) != regionTriangles.end() )
 326                 regionTriangles.insert
 327                 (
 328                     std::make_pair
 329                     (
 330                         regionI,
 331                         DynList<label>()
 332                     )
 333                 );
 334             if( !otherLabels.found(regionI) )
 335                 otherLabels.insert(regionI, labelHashSet());
 336
 337             regionTriangles[regionI].append(triI);
 338
 339             forAll(surface_[triI], tpI)
 340             {
 341                 const label pI = surface_[triI][tpI];
 342
 343                 if( pointI == pI )
 344                     continue;
 345
 346                 otherLabels[regionI].insert(pI);
 347                 normals[regionI] += surface_[triI].normal(points);
 348             }
 349         }
 350
 351         forAllConstIter(Map<labelHashSet>, otherLabels, it)
 352         {
 353             const labelHashSet& currLabels = it();
 354
 355             if( otherLabels[it.key()].size() > 5 )
 356                 continue;
 357
 358             labelHashSet additionalPoints;
 359             forAllConstIter(labelHashSet, currLabels, lit)
 360             {
 361                 const label neiPointI = lit.key();
 362                 const constRow pTriangles = pointTriangles[neiPointI];
 363
 364                 forAll(pTriangles, ptI)
 365                 {
 366                     const labelledTri& nTri = surface_[pTriangles[ptI]];
 367
 368                     if( nTri.region() != it.key() )
 369                         continue;
 370
 371                     forAll(nTri, pI)
 372                         additionalPoints.insert(nTri[pI]);
 373                 }
 374             }
 375
 376             forAllConstIter(labelHashSet, additionalPoints, aIter)
 377                 otherLabels[it.key()].insert(aIter.key());
 378         }
 379
 380         forAllIter(Map<vector>, normals, nit)
 381             nit() /= (Foam::mag(nit()) + VSMALL);
 382
 383         forAllConstIter(Map<labelHashSet>, otherLabels, it)
 384         {
 385             const labelHashSet& labels = it();
 386
 387             const label regionI = it.key();
 388
 389             //- store the patches
 390             pointPatches[pointI].append(regionI);
 391
 392             //- find the position of point in the patch triangles
 393             const DynList<label>& rTriangles = regionTriangles[regionI];
 394             forAll(rTriangles, i)
 395             {
 396                 const label tI = rTriangles[i];
 397
 398                 label pos(-1);
 399                 forAll(surface_[tI], j)
 400                     if( surface_[tI][j] == pointI )
 401                     {
 402                         pos = j;
 403                         break;
 404                     }
 405
 406                 patchPositions_(tI, pos) = gaussianCurvature_[pointI].size();
 407             }
 408
 409             //- store point coordinates
 410             DynList<point> op;
 411
 412             forAllConstIter(labelHashSet, labels, lit)
 413                 op.append(points[lit.key()]);
 414
 415             //- fit the quadric patch to the surface
 416             quadricFitting qfit(points[pointI], normals[it.key()], op);
 417
 418             # ifdef DEBUGCurvatureEstimator
 419             Info << "Point " << pointI << " in patch " << regionI
 420                 << " normal " << normals[it.key()]
 421                 << " evalution points " << labels
 422                 << " has max curvature " << qfit.maxCurvature()
 423                 << " and min curvature " << qfit.minCurvature() << endl;
 424             OFstream file
 425             (
 426                 "point_" +
 427                 help::scalarToText(pointI) +
 428                 "_region_" +
 429                 help::scalarToText(regionI) +
 430                 "_triangles.vtk"
 431             );
 432             writeSurfaceToVTK(file, surface_, rTriangles, labels);
 433             # endif
 434
 435             //- store curvatures
 436             gaussianCurvature_[pointI].append(qfit.gaussianCurvature());
 437             meanCurvature_[pointI].append(qfit.meanCurvature());
 438             maxCurvature_[pointI].append(qfit.maxCurvature());
 439             minCurvature_[pointI].append(qfit.minCurvature());
 440             maxCurvatureVector_[pointI].append(qfit.maxCurvatureVector());
 441             minCurvatureVector_[pointI].append(qfit.minCurvatureVector());
 442         }
 443     }
 444
 445     //- smooth curvatures using weighted Laplace
 446     List<DynList<scalar, 1> > smoothMinCurv(points.size());
 447     List<DynList<scalar, 1> > smoothMaxCurv(points.size());
 448
 449     # ifdef USE_OMP
 450     # pragma omp parallel
 451     # endif
 452     {
 453         for(label iteration=0;iteration<2;++iteration)
 454         {
 455             # ifdef USE_OMP
 456             # pragma omp for schedule(static, 1)
 457             # endif
 458             forAll(pointTriangles, pointI)
 459             {
 460                 const constRow pTriangles = pointTriangles[pointI];
 461
 462                 //- find neighbouring points for each patch
 463                 Map<DynList<label> > patchNeiPoints;
 464                 forAll(pointPatches[pointI], ppI)
 465                     patchNeiPoints.insert
 466                     (
 467                         pointPatches[pointI][ppI],
 468                         DynList<label>()
 469                     );
 470
 471                 forAll(pTriangles, ptI)
 472                 {
 473                     const labelledTri& tri = surface_[pTriangles[ptI]];
 474
 475                     if( !patchNeiPoints.found(tri.region()) )
 476                         patchNeiPoints.insert(tri.region(), DynList<label>());
 477
 478                     forAll(tri, pI)
 479                     {
 480                         const label neiPointI = tri[pI];
 481
 482                         if( neiPointI == pointI )
 483                             continue;
 484
 485                         patchNeiPoints[tri.region()].appendIfNotIn(neiPointI);
 486                     }
 487                 }
 488
 489                 smoothMinCurv[pointI].setSize(minCurvature_[pointI].size());
 490                 smoothMaxCurv[pointI].setSize(maxCurvature_[pointI].size());
 491
 492                 //- update min curvature for all point patches
 493                 forAll(minCurvature_[pointI], patchI)
 494                 {
 495                     const label cPatch = pointPatches[pointI][patchI];
 496
 497                     scalar minCurv(0.0);
 498                     scalar maxCurv(0.0);
 499
 500                     const DynList<label>& neiPoints = patchNeiPoints[cPatch];
 501
 502                     if( neiPoints.size() == 0 )
 503                     {
 504                         smoothMinCurv[pointI][patchI] =
 505                             minCurvature_[pointI][patchI];
 506
 507                         smoothMaxCurv[pointI][patchI] =
 508                             maxCurvature_[pointI][patchI];
 509                     }
 510
 511                     forAll(neiPoints, i)
 512                     {
 513                         const label neiPointI = neiPoints[i];
 514
 515                         const label pos =
 516                             pointPatches[neiPointI].containsAtPosition(cPatch);
 517
 518                         minCurv += minCurvature_[neiPointI][pos];
 519                         maxCurv += maxCurvature_[neiPointI][pos];
 520                     }
 521
 522                     minCurv /= neiPoints.size();
 523                     maxCurv /= neiPoints.size();
 524
 525                     //- store the value
 526                     smoothMinCurv[pointI][patchI] =
 527                         0.5 * (minCurv + minCurvature_[pointI][patchI]);
 528
 529                     smoothMaxCurv[pointI][patchI] =
 530                         0.5 * (maxCurv + maxCurvature_[pointI][patchI]);
 531                 }
 532             }
 533
 534             # ifdef USE_OMP
 535             # pragma omp barrier
 536
 537             # pragma omp for schedule(static, 1)
 538             # endif
 539             forAll(minCurvature_, pointI)
 540             {
 541                 forAll(minCurvature_[pointI], i)
 542                 {
 543                     minCurvature_[pointI][i] = smoothMinCurv[pointI][i];
 544                     maxCurvature_[pointI][i] = smoothMaxCurv[pointI][i];
 545                 }
 546             }
 547         }
 548
 549         # ifdef USE_OMP
 550         # pragma omp barrier
 551         # endif
 552
 553         //- update Gaussian and mean curvatures
 554         # ifdef USE_OMP
 555         # pragma omp for schedule(static, 1)
 556         # endif
 557         forAll(minCurvature_, pointI)
 558         {
 559             const DynList<scalar, 1>& minCurv = minCurvature_[pointI];
 560             const DynList<scalar, 1>& maxCurv = maxCurvature_[pointI];
 561
 562             forAll(minCurv, i)
 563             {
 564                 gaussianCurvature_[pointI][i] = minCurv[i] * maxCurv[i];
 565                 meanCurvature_[pointI][i] = 0.5 * (minCurv[i] + maxCurv[i]);
 566             }
 567         }
 568     }
 569
 570     # ifdef DEBUGCurvatureEstimator
 571     word name = "surfaceMeanCurv";
 572     writeSurfaceToVTK(name, surface_, meanCurvature_);
 573     writeSurfaceToVTK("surfaceGaussianCurv", surface_, gaussianCurvature_);
 574     writeSurfaceToVTK("surfaceMaxCurv", surface_, maxCurvature_);
 575     writeSurfaceToVTK("surfaceMinCurv", surface_, minCurvature_);
 576     writeSurfaceToVTK("surfaceMaxCurvVec", surface_, maxCurvatureVector_);
 577     writeSurfaceToVTK("surfaceMinCurvVec", surface_, minCurvatureVector_);
 578     # endif
 579 }
 580 /*
 581 void triSurfaceCurvatureEstimator::calculateGaussianCurvature()
 582 {
 583     gaussianCurvature_.setSize(surface_.patches().size());
 584
 585     List<std::map<label, scalar> > magAreas(gaussianCurvature_.size());
 586
 587     forAll(surface_, triI)
 588     {
 589         const labelledTri& tri = surface_[triI];
 590         const label region = tri.region();
 591
 592         forAll(tri, pI)
 593         {
 594             gaussianCurvature_[region][tri[pI]] = 2 * M_PI;
 595             magAreas[region][tri[pI]] = 0.0;
 596         }
 597     }
 598
 599     const pointField& points = surface_.points();
 600     const labelListList& pointTriangles = surface_.pointFaces();
 601
 602     forAll(surface_, triI)
 603     {
 604         const labelledTri& tri = surface_[triI];
 605         const label region = tri.region();
 606
 607         const scalar A = mag(tri.normal(points)) + VSMALL;
 608
 609         const edgeList edges = tri.edges();
 610         vectorField ev(3);
 611         scalarField dv(3);
 612         forAll(edges, eI)
 613         {
 614             ev[eI] = edges[eI].vec(points);
 615             dv[eI] = mag(ev[eI]);
 616             if( dv[eI] > VSMALL )
 617                 ev[eI] /= dv[eI];
 618         }
 619
 620         forAll(tri, pI)
 621         {
 622             scalar cs = ev[(pI+2)%3] & ev[pI];
 623             cs = Foam::min(1.0, cs);
 624             cs = Foam::max(1.0, cs);
 625
 626             gaussianCurvature_[region][tri[pI]] -= Foam::acos(-1.0 * cs);
 627             magAreas[region][tri[pI]] += A;
 628         }
 629     }
 630
 631     //- calculate tge gaussian curvature for each vertex
 632     forAll(gaussianCurvature_, patchI)
 633     {
 634         std::map<label, scalar>& gc = gaussianCurvature_[patchI];
 635         std::map<label, scalar>& magSf = magAreas[patchI];
 636
 637         for
 638         (
 639             std::map<label, scalar>::iterator it = gc.begin();
 640             it!=gc.end();
 641             ++it
 642         )
 643             it->second = 3.0 * it->second / magSf[it->first];
 644     }
 645
 646     //- smooth the curvature variation
 647     const edgeList& edges = surface_.edges();
 648     const labelListList& pointEdges = surface_.pointEdges();
 649     forAll(meanCurvature_, patchI)
 650     {
 651         std::map<label, scalar>& gc = meanCurvature_[patchI];
 652
 653         std::map<label, scalar> newValues;
 654
 655         for(std::map<label, scalar>::iterator it=gc.begin();it!=gc.end();++it)
 656         {
 657             const labelList& pEdges = pointEdges[it->first];
 658
 659             scalar maxCurv(-VSMALL), minCurv(VSMALL);
 660             label nNei(0);
 661
 662             forAll(pEdges, peI)
 663             {
 664                 const label ovI = edges[pEdges[peI]].otherVertex(it->first);
 665
 666                 if( gc.find(ovI) != gc.end() )
 667                 {
 668                     maxCurv = Foam::max(maxCurv, gc[ovI]);
 669                     minCurv = Foam::min(minCurv, gc[ovI]);
 670                     ++nNei;
 671                 }
 672             }
 673
 674             if( nNei != 0 )
 675             {
 676                 if( gc[it->first] > maxCurv )
 677                 {
 678                     newValues.insert
 679                     (
 680                         std::make_pair(it->first, 0.5*(maxCurv+gc[it->first]))
 681                     );
 682                 }
 683                 else if( gc[it->first] >= minCurv )
 684                 {
 685                     newValues.insert(*it);
 686                 }
 687                 else
 688                 {
 689                     newValues.insert
 690                     (
 691                         std::make_pair(it->first, 0.5*(minCurv+gc[it->first]))
 692                     );
 693                 }
 694             }
 695         }
 696
 697         gc = newValues;
 698     }
 699 }
 700
 701 void triSurfaceCurvatureEstimator::calculateMeanCurvature()
 702 {
 703     const pointField& points = surface_.points();
 704
 705     meanCurvature_.setSize(surface_.patches().size());
 706
 707     List<Map<label> > nNeighbours(meanCurvature_.size());
 708
 709     const edgeList& edges = surface_.edges();
 710     const labelListList& edgeFaces = surface_.edgeFaces();
 711
 712     //- calculate edge angles
 713     forAll(edges, eI)
 714     {
 715         const labelList& eFaces = edgeFaces[eI];
 716
 717         if( eFaces.size() != 2 )
 718             continue;
 719         if( surface_[eFaces[0]].region() != surface_[eFaces[1]].region() )
 720             continue;
 721
 722         const label region = surface_[eFaces[0]].region();
 723         std::map<label, scalar>& mc = meanCurvature_[region];
 724         Map<label>& nn = nNeighbours[region];
 725
 726         //- calculate normal vectors
 727         vector n0 = surface_[eFaces[0]].normal(points);
 728         const scalar A0 = mag(n0);
 729         if( A0 > VSMALL )
 730             n0 /= A0;
 731
 732         vector n1 = surface_[eFaces[1]].normal(points);
 733         const scalar A1 = mag(n1);
 734         if( A1 > VSMALL )
 735             n1 /= A1;
 736
 737         //- calculate the cosine of the angle between the two vectors
 738         scalar cs = n0 & n1;
 739         cs = Foam::min(1.0, cs);
 740         cs = Foam::max(-1.0, cs);
 741
 742         //- calculate and normalize the edge vector
 743         vector ev = edges[eI].vec(points);
 744         const scalar length = mag(ev);
 745         if( length > VSMALL )
 746             ev /= length;
 747
 748         //- calculate cross product of the normal vectors
 749         const vector t = n0 ^ n1;
 750
 751         //- calculate the dot product between the edge vector and t
 752         scalar sn = t & ev;
 753         sn = Foam::min(1.0, sn);
 754         sn = Foam::max(-1.0, sn);
 755
 756         //- calculate mean curvature over the edge
 757         scalar Hf(0.0);
 758
 759         if( mag(sn) > VSMALL || mag(cs) > VSMALL )
 760         {
 761             //- calculate the angle
 762             const scalar angle = Foam::atan2(sn, cs);
 763             Hf += angle * length;
 764         }
 765
 766         if( A0 > VSMALL || A1 > VSMALL )
 767         {
 768             Hf /= (A0 + A1 + VSMALL);
 769             Hf *= 3.0;
 770         }
 771
 772         //- store the data
 773         const edge& e = edges[eI];
 774         if( mc.find(e.start()) == mc.end() )
 775         {
 776             mc.insert(std::make_pair(e.start(), 0.0));
 777             nn.insert(e.start(), 0);
 778         }
 779         mc[e.start()] += Hf;
 780         nn[e.start()] += 1;
 781
 782         if( mc.find(e.end()) == mc.end() )
 783         {
 784             mc.insert(std::make_pair(e.end(), 0.0));
 785             nn.insert(e.end(), 0);
 786         }
 787         mc[e.end()] += Hf;
 788         nn[e.end()] += 1;
 789     }
 790
 791     //- calculate the curvature for each vertex
 792     forAll(meanCurvature_, patchI)
 793     {
 794         std::map<label, scalar>& mc = meanCurvature_[patchI];
 795         Map<label>& nn = nNeighbours[patchI];
 796
 797         forAllIter(Map<label>, nn, it)
 798         {
 799             if( it() == 0 )
 800             {
 801                 mc[it.key()] = 0.0;
 802                 continue;
 803             }
 804
 805             mc[it.key()] = 0.5 * mc[it.key()] / it();
 806         }
 807     }
 808
 809     //- smooth the curvature variation
 810     const labelListList& pointEdges = surface_.pointEdges();
 811     forAll(meanCurvature_, patchI)
 812     {
 813         std::map<label, scalar>& mc = meanCurvature_[patchI];
 814
 815         std::map<label, scalar> newValues;
 816
 817         for(std::map<label, scalar>::iterator it=mc.begin();it!=mc.end();++it)
 818         {
 819             const labelList& pEdges = pointEdges[it->first];
 820
 821             scalar maxCurv(-VSMALL), minCurv(VSMALL);
 822             label nNei(0);
 823
 824             forAll(pEdges, peI)
 825             {
 826                 const label ovI = edges[pEdges[peI]].otherVertex(it->first);
 827
 828                 if( mc.find(ovI) != mc.end() )
 829                 {
 830                     maxCurv = Foam::max(maxCurv, mc[ovI]);
 831                     minCurv = Foam::min(minCurv, mc[ovI]);
 832                     ++nNei;
 833                 }
 834             }
 835
 836             if( nNei != 0 )
 837             {
 838                 if( mc[it->first] > maxCurv )
 839                 {
 840                     newValues.insert
 841                     (
 842                         std::make_pair(it->first, 0.5*(maxCurv+mc[it->first]))
 843                     );
 844                 }
 845                 else if( mc[it->first] >= minCurv )
 846                 {
 847                     newValues.insert(*it);
 848                 }
 849                 else
 850                 {
 851                     newValues.insert
 852                     (
 853                         std::make_pair(it->first, 0.5*(minCurv+mc[it->first]))
 854                     );
 855                 }
 856             }
 857         }
 858
 859         mc = newValues;
 860     }
 861 }
 862
 863 void triSurfaceCurvatureEstimator::calculateMinAndMaxCurvature()
 864 {
 865     maxCurvature_.setSize(meanCurvature_.size());
 866     minCurvature_.setSize(meanCurvature_.size());
 867
 868     forAll(maxCurvature_, patchI)
 869     {
 870         std::map<label, scalar>& mc = meanCurvature_[patchI];
 871         std::map<label, scalar>& gc = gaussianCurvature_[patchI];
 872         for
 873         (
 874             std::map<label, scalar>::const_iterator it=mc.begin();
 875             it!=mc.end();
 876             ++it
 877         )
 878         {
 879             const label pI = it->first;
 880             scalar disc = sqr(mc[pI]) - gc[pI];
 881             if( disc <  0 )
 882                 disc = sqr(mc[pI]) + gc[pI];
 883
 884             const scalar sqrtdisc = sqrt(disc);
 885             maxCurvature_[patchI][pI] = mc[pI] + sqrtdisc;
 886             minCurvature_[patchI][pI] = mc[pI] - sqrtdisc;
 887         }
 888     }
 889 }
 890 */
 891
 892 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 893
 894 } // End namespace Foam
 895
 896 // ************************************************************************* //