src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "extendedFeatureEdgeMesh.H"
  27 #include "triSurface.H"
  28 #include "Random.H"
  29 #include "Time.H"
  30 #include "meshTools.H"
  31 #include "linePointRef.H"
  32 #include "ListListOps.H"
  33 #include "OFstream.H"
  34 #include "IFstream.H"
  35 #include "unitConversion.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(Foam::extendedFeatureEdgeMesh, 0);
  40
  41 Foam::scalar Foam::extendedFeatureEdgeMesh::cosNormalAngleTol_ =
  42     Foam::cos(degToRad(0.1));
  43
  44
  45 Foam::label Foam::extendedFeatureEdgeMesh::convexStart_ = 0;
  46
  47
  48 Foam::label Foam::extendedFeatureEdgeMesh::externalStart_ = 0;
  49
  50
  51 Foam::label Foam::extendedFeatureEdgeMesh::nPointTypes = 4;
  52
  53
  54 Foam::label Foam::extendedFeatureEdgeMesh::nEdgeTypes = 5;
  55
  56
  57 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  58
  59 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
  60 :
  61     regIOobject(io),
  62     edgeMesh(pointField(0), edgeList(0)),
  63     concaveStart_(0),
  64     mixedStart_(0),
  65     nonFeatureStart_(0),
  66     internalStart_(0),
  67     flatStart_(0),
  68     openStart_(0),
  69     multipleStart_(0),
  70     normals_(0),
  71     edgeDirections_(0),
  72     edgeNormals_(0),
  73     featurePointNormals_(0),
  74     regionEdges_(0),
  75     edgeTree_(),
  76     edgeTreesByType_()
  77 {
  78     if
  79     (
  80         io.readOpt() == IOobject::MUST_READ
  81      || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
  82      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
  83     )
  84     {
  85         if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
  86         {
  87             WarningIn
  88             (
  89                 "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
  90                 "(const IOobject&)"
  91             )   << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
  92                 << " does not support automatic rereading."
  93                 << endl;
  94         }
  95
  96         Istream& is = readStream(typeName);
  97
  98         is  >> *this
  99             >> concaveStart_
 100             >> mixedStart_
 101             >> nonFeatureStart_
 102             >> internalStart_
 103             >> flatStart_
 104             >> openStart_
 105             >> multipleStart_
 106             >> normals_
 107             >> edgeNormals_
 108             >> featurePointNormals_
 109             >> regionEdges_;
 110
 111         close();
 112
 113         {
 114             // Calculate edgeDirections
 115
 116             const edgeList& eds(edges());
 117
 118             const pointField& pts(points());
 119
 120             edgeDirections_.setSize(eds.size());
 121
 122             forAll(eds, eI)
 123             {
 124                 edgeDirections_[eI] = eds[eI].vec(pts);
 125             }
 126
 127             edgeDirections_ /= mag(edgeDirections_);
 128         }
 129     }
 130
 131     if (debug)
 132     {
 133         Pout<< "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh :"
 134             << " constructed from IOobject :"
 135             << " points:" << points().size()
 136             << " edges:" << edges().size()
 137             << endl;
 138     }
 139 }
 140
 141
 142 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 143 (
 144     const IOobject& io,
 145     const extendedFeatureEdgeMesh& fem
 146 )
 147 :
 148     regIOobject(io),
 149     edgeMesh(fem),
 150     concaveStart_(fem.concaveStart()),
 151     mixedStart_(fem.mixedStart()),
 152     nonFeatureStart_(fem.nonFeatureStart()),
 153     internalStart_(fem.internalStart()),
 154     flatStart_(fem.flatStart()),
 155     openStart_(fem.openStart()),
 156     multipleStart_(fem.multipleStart()),
 157     normals_(fem.normals()),
 158     edgeDirections_(fem.edgeDirections()),
 159     edgeNormals_(fem.edgeNormals()),
 160     featurePointNormals_(fem.featurePointNormals()),
 161     regionEdges_(fem.regionEdges()),
 162     edgeTree_(),
 163     edgeTreesByType_()
 164 {}
 165
 166
 167 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 168 (
 169     const IOobject& io,
 170     const Xfer<pointField>& pointLst,
 171     const Xfer<edgeList>& edgeLst
 172 )
 173 :
 174     regIOobject(io),
 175     edgeMesh(pointLst, edgeLst),
 176     concaveStart_(0),
 177     mixedStart_(0),
 178     nonFeatureStart_(0),
 179     internalStart_(0),
 180     flatStart_(0),
 181     openStart_(0),
 182     multipleStart_(0),
 183     normals_(0),
 184     edgeDirections_(0),
 185     edgeNormals_(0),
 186     featurePointNormals_(0),
 187     regionEdges_(0),
 188     edgeTree_(),
 189     edgeTreesByType_()
 190 {}
 191
 192
 193 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 194 (
 195     const surfaceFeatures& sFeat,
 196     const objectRegistry& obr,
 197     const fileName& sFeatFileName
 198 )
 199 :
 200     regIOobject
 201     (
 202         IOobject
 203         (
 204             sFeatFileName,
 205             obr.time().constant(),
 206             "extendedFeatureEdgeMesh",
 207             obr,
 208             IOobject::NO_READ,
 209             IOobject::NO_WRITE
 210         )
 211     ),
 212     edgeMesh(pointField(0), edgeList(0)),
 213     concaveStart_(-1),
 214     mixedStart_(-1),
 215     nonFeatureStart_(-1),
 216     internalStart_(-1),
 217     flatStart_(-1),
 218     openStart_(-1),
 219     multipleStart_(-1),
 220     normals_(0),
 221     edgeDirections_(0),
 222     edgeNormals_(0),
 223     featurePointNormals_(0),
 224     regionEdges_(0),
 225     edgeTree_(),
 226     edgeTreesByType_()
 227 {
 228     // Extract and reorder the data from surfaceFeatures
 229
 230     // References to the surfaceFeatures data
 231     const triSurface& surf(sFeat.surface());
 232     const pointField& sFeatLocalPts(surf.localPoints());
 233     const edgeList& sFeatEds(surf.edges());
 234
 235     // Filling the extendedFeatureEdgeMesh with the raw geometrical data.
 236
 237     label nFeatEds = sFeat.featureEdges().size();
 238
 239     DynamicList<point> tmpPts;
 240     edgeList eds(nFeatEds);
 241     DynamicList<vector> norms;
 242     vectorField edgeDirections(nFeatEds);
 243     labelListList edgeNormals(nFeatEds);
 244     DynamicList<label> regionEdges;
 245
 246     // Mapping between old and new indices, there is entry in the map for each
 247     // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
 248     // corresponds to the index
 249     labelList pointMap(sFeatLocalPts.size(), -1);
 250
 251     // Noting when the normal of a face has been used so not to duplicate
 252     labelList faceMap(surf.size(), -1);
 253
 254     // Collecting the status of edge for subsequent sorting
 255     List<edgeStatus> edStatus(nFeatEds, NONE);
 256
 257     forAll(sFeat.featurePoints(), i)
 258     {
 259         label sFPI = sFeat.featurePoints()[i];
 260
 261         tmpPts.append(sFeatLocalPts[sFPI]);
 262
 263         pointMap[sFPI] = tmpPts.size() - 1;
 264     }
 265
 266     // All feature points have been added
 267     nonFeatureStart_ = tmpPts.size();
 268
 269     forAll(sFeat.featureEdges(), i)
 270     {
 271         label sFEI = sFeat.featureEdges()[i];
 272
 273         const edge& fE(sFeatEds[sFEI]);
 274
 275         // Check to see if the points have been already used
 276         if (pointMap[fE.start()] == -1)
 277         {
 278              tmpPts.append(sFeatLocalPts[fE.start()]);
 279
 280              pointMap[fE.start()] = tmpPts.size() - 1;
 281         }
 282
 283         eds[i].start() = pointMap[fE.start()];
 284
 285         if (pointMap[fE.end()] == -1)
 286         {
 287             tmpPts.append(sFeatLocalPts[fE.end()]);
 288
 289             pointMap[fE.end()] = tmpPts.size() - 1;
 290         }
 291
 292         eds[i].end() = pointMap[fE.end()];
 293
 294         // Pick up the faces adjacent to the feature edge
 295         const labelList& eFaces = surf.edgeFaces()[sFEI];
 296
 297         edgeNormals[i].setSize(eFaces.size());
 298
 299         forAll(eFaces, j)
 300         {
 301             label eFI = eFaces[j];
 302
 303             // Check to see if the points have been already used
 304             if (faceMap[eFI] == -1)
 305             {
 306                 norms.append(surf.faceNormals()[eFI]);
 307
 308                 faceMap[eFI] = norms.size() - 1;
 309             }
 310
 311             edgeNormals[i][j] = faceMap[eFI];
 312         }
 313
 314         vector fC0tofC1(vector::zero);
 315
 316         if (eFaces.size() == 2)
 317         {
 318             fC0tofC1 =
 319                 surf[eFaces[1]].centre(surf.points())
 320               - surf[eFaces[0]].centre(surf.points());
 321         }
 322
 323         edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
 324
 325         edgeDirections[i] = fE.vec(sFeatLocalPts);
 326
 327         if (i < sFeat.nRegionEdges())
 328         {
 329             regionEdges.append(i);
 330         }
 331     }
 332
 333     // Reorder the edges by classification
 334
 335     List<DynamicList<label> > allEds(nEdgeTypes);
 336
 337     DynamicList<label>& externalEds(allEds[0]);
 338     DynamicList<label>& internalEds(allEds[1]);
 339     DynamicList<label>& flatEds(allEds[2]);
 340     DynamicList<label>& openEds(allEds[3]);
 341     DynamicList<label>& multipleEds(allEds[4]);
 342
 343     forAll(eds, i)
 344     {
 345         edgeStatus eStat = edStatus[i];
 346
 347         if (eStat == EXTERNAL)
 348         {
 349             externalEds.append(i);
 350         }
 351         else if (eStat == INTERNAL)
 352         {
 353             internalEds.append(i);
 354         }
 355         else if (eStat == FLAT)
 356         {
 357             flatEds.append(i);
 358         }
 359         else if (eStat == OPEN)
 360         {
 361             openEds.append(i);
 362         }
 363         else if (eStat == MULTIPLE)
 364         {
 365             multipleEds.append(i);
 366         }
 367         else if (eStat == NONE)
 368         {
 369             FatalErrorIn
 370             (
 371                 "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
 372                 "("
 373                 "    const surfaceFeatures& sFeat,"
 374                 "    const objectRegistry& obr,"
 375                 "    const fileName& sFeatFileName"
 376                 ")"
 377             )
 378                 << nl << "classifyEdge returned NONE on edge "
 379                 << eds[i]
 380                 << ". There is a problem with definition of this edge."
 381                 << nl << abort(FatalError);
 382         }
 383     }
 384
 385     internalStart_ = externalEds.size();
 386     flatStart_ = internalStart_ + internalEds.size();
 387     openStart_ = flatStart_ + flatEds.size();
 388     multipleStart_ = openStart_ + openEds.size();
 389
 390     labelList edMap
 391     (
 392         ListListOps::combine<labelList>
 393         (
 394             allEds,
 395             accessOp<labelList>()
 396         )
 397     );
 398
 399     edMap = invert(edMap.size(), edMap);
 400
 401     inplaceReorder(edMap, eds);
 402     inplaceReorder(edMap, edStatus);
 403     inplaceReorder(edMap, edgeDirections);
 404     inplaceReorder(edMap, edgeNormals);
 405     inplaceRenumber(edMap, regionEdges);
 406
 407     pointField pts(tmpPts);
 408
 409     // Initialise the edgeMesh
 410     edgeMesh::operator=(edgeMesh(pts, eds));
 411
 412     // Initialise sorted edge related data
 413     edgeDirections_ = edgeDirections/mag(edgeDirections);
 414     edgeNormals_ = edgeNormals;
 415     regionEdges_ = regionEdges;
 416
 417     // Normals are all now found and indirectly addressed, can also be stored
 418     normals_ = vectorField(norms);
 419
 420     // Reorder the feature points by classification
 421
 422     List<DynamicList<label> > allPts(3);
 423
 424     DynamicList<label>& convexPts(allPts[0]);
 425     DynamicList<label>& concavePts(allPts[1]);
 426     DynamicList<label>& mixedPts(allPts[2]);
 427
 428     for (label i = 0; i < nonFeatureStart_; i++)
 429     {
 430         pointStatus ptStatus = classifyFeaturePoint(i);
 431
 432         if (ptStatus == CONVEX)
 433         {
 434             convexPts.append(i);
 435         }
 436         else if (ptStatus == CONCAVE)
 437         {
 438             concavePts.append(i);
 439         }
 440         else if (ptStatus == MIXED)
 441         {
 442             mixedPts.append(i);
 443         }
 444         else if (ptStatus == NONFEATURE)
 445         {
 446             FatalErrorIn
 447             (
 448                 "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
 449                 "("
 450                 "    const surfaceFeatures& sFeat,"
 451                 "    const objectRegistry& obr,"
 452                 "    const fileName& sFeatFileName"
 453                 ")"
 454             )
 455                 << nl << "classifyFeaturePoint returned NONFEATURE on point at "
 456                 << points()[i]
 457                 << ". There is a problem with definition of this feature point."
 458                 << nl << abort(FatalError);
 459         }
 460     }
 461
 462     concaveStart_ = convexPts.size();
 463     mixedStart_ = concaveStart_ + concavePts.size();
 464
 465     labelList ftPtMap
 466     (
 467         ListListOps::combine<labelList>
 468         (
 469             allPts,
 470             accessOp<labelList>()
 471         )
 472     );
 473
 474     ftPtMap = invert(ftPtMap.size(), ftPtMap);
 475
 476     // Creating the ptMap from the ftPtMap with identity values up to the size
 477     // of pts to create an oldToNew map for inplaceReorder
 478
 479     labelList ptMap(identity(pts.size()));
 480
 481     forAll(ftPtMap, i)
 482     {
 483         ptMap[i] = ftPtMap[i];
 484     }
 485
 486     inplaceReorder(ptMap, pts);
 487
 488     forAll(eds, i)
 489     {
 490         inplaceRenumber(ptMap, eds[i]);
 491     }
 492
 493     // Reinitialise the edgeMesh with sorted feature points and
 494     // renumbered edges
 495     edgeMesh::operator=(edgeMesh(pts, eds));
 496
 497     // Generate the featurePointNormals
 498
 499     labelListList featurePointNormals(nonFeatureStart_);
 500
 501     for (label i = 0; i < nonFeatureStart_; i++)
 502     {
 503         DynamicList<label> tmpFtPtNorms;
 504
 505         const labelList& ptEds = pointEdges()[i];
 506
 507         forAll(ptEds, j)
 508         {
 509             const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
 510
 511             forAll(ptEdNorms, k)
 512             {
 513                 if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
 514                 {
 515                     bool addNormal = true;
 516
 517                     // Check that the normal direction is unique at this feature
 518                     forAll(tmpFtPtNorms, q)
 519                     {
 520                         if
 521                         (
 522                             (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
 523                           > cosNormalAngleTol_
 524                         )
 525                         {
 526                             // Parallel to an existing normal, do not add
 527                             addNormal = false;
 528
 529                             break;
 530                         }
 531                     }
 532
 533                     if (addNormal)
 534                     {
 535                         tmpFtPtNorms.append(ptEdNorms[k]);
 536                     }
 537                 }
 538             }
 539         }
 540
 541         featurePointNormals[i] = tmpFtPtNorms;
 542     }
 543
 544     featurePointNormals_ = featurePointNormals;
 545 }
 546
 547
 548 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 549 (
 550     const IOobject& io,
 551     const pointField& pts,
 552     const edgeList& eds,
 553     label concaveStart,
 554     label mixedStart,
 555     label nonFeatureStart,
 556     label internalStart,
 557     label flatStart,
 558     label openStart,
 559     label multipleStart,
 560     const vectorField& normals,
 561     const vectorField& edgeDirections,
 562     const labelListList& edgeNormals,
 563     const labelListList& featurePointNormals,
 564     const labelList& regionEdges
 565 )
 566 :
 567     regIOobject(io),
 568     edgeMesh(pts, eds),
 569     concaveStart_(concaveStart),
 570     mixedStart_(mixedStart),
 571     nonFeatureStart_(nonFeatureStart),
 572     internalStart_(internalStart),
 573     flatStart_(flatStart),
 574     openStart_(openStart),
 575     multipleStart_(multipleStart),
 576     normals_(normals),
 577     edgeDirections_(edgeDirections),
 578     edgeNormals_(edgeNormals),
 579     featurePointNormals_(featurePointNormals),
 580     regionEdges_(regionEdges),
 581     edgeTree_(),
 582     edgeTreesByType_()
 583 {}
 584
 585
 586 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 587
 588 Foam::extendedFeatureEdgeMesh::~extendedFeatureEdgeMesh()
 589 {}
 590
 591
 592 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 593
 594 Foam::extendedFeatureEdgeMesh::pointStatus
 595 Foam::extendedFeatureEdgeMesh::classifyFeaturePoint
 596 (
 597     label ptI
 598 ) const
 599 {
 600     labelList ptEds(pointEdges()[ptI]);
 601
 602     label nPtEds = ptEds.size();
 603     label nExternal = 0;
 604     label nInternal = 0;
 605
 606     if (nPtEds == 0)
 607     {
 608         // There are no edges attached to the point, this is a problem
 609         return NONFEATURE;
 610     }
 611
 612     forAll(ptEds, i)
 613     {
 614         edgeStatus edStat = getEdgeStatus(ptEds[i]);
 615
 616         if (edStat == EXTERNAL)
 617         {
 618             nExternal++;
 619         }
 620         else if (edStat == INTERNAL)
 621         {
 622             nInternal++;
 623         }
 624     }
 625
 626     if (nExternal == nPtEds)
 627     {
 628         return CONVEX;
 629     }
 630     else if (nInternal == nPtEds)
 631     {
 632         return CONCAVE;
 633     }
 634     else
 635     {
 636         return MIXED;
 637     }
 638 }
 639
 640
 641 Foam::extendedFeatureEdgeMesh::edgeStatus
 642 Foam::extendedFeatureEdgeMesh::classifyEdge
 643 (
 644     const List<vector>& norms,
 645     const labelList& edNorms,
 646     const vector& fC0tofC1
 647 ) const
 648 {
 649     label nEdNorms = edNorms.size();
 650
 651     if (nEdNorms == 1)
 652     {
 653         return OPEN;
 654     }
 655     else if (nEdNorms == 2)
 656     {
 657         const vector n0(norms[edNorms[0]]);
 658         const vector n1(norms[edNorms[1]]);
 659
 660         if ((n0 & n1) > cosNormalAngleTol_)
 661         {
 662             return FLAT;
 663         }
 664         else if ((fC0tofC1 & n0) > 0.0)
 665         {
 666             return INTERNAL;
 667         }
 668         else
 669         {
 670             return EXTERNAL;
 671         }
 672     }
 673     else if (nEdNorms > 2)
 674     {
 675         return MULTIPLE;
 676     }
 677     else
 678     {
 679         // There is a problem - the edge has no normals
 680         return NONE;
 681     }
 682 }
 683
 684
 685 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 686
 687 void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
 688 (
 689     const point& sample,
 690     scalar searchDistSqr,
 691     pointIndexHit& info
 692 ) const
 693 {
 694     info = edgeTree().findNearest
 695     (
 696         sample,
 697         searchDistSqr
 698     );
 699 }
 700
 701
 702 void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
 703 (
 704     const pointField& samples,
 705     const scalarField& searchDistSqr,
 706     List<pointIndexHit>& info
 707 ) const
 708 {
 709     info.setSize(samples.size());
 710
 711     forAll(samples, i)
 712     {
 713         nearestFeatureEdge
 714         (
 715             samples[i],
 716             searchDistSqr[i],
 717             info[i]
 718         );
 719     }
 720 }
 721
 722
 723 void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType
 724 (
 725     const point& sample,
 726     const scalarField& searchDistSqr,
 727     List<pointIndexHit>& info
 728 ) const
 729 {
 730     const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
 731
 732     info.setSize(edgeTrees.size());
 733
 734     labelList sliceStarts(edgeTrees.size());
 735
 736     sliceStarts[0] = externalStart_;
 737     sliceStarts[1] = internalStart_;
 738     sliceStarts[2] = flatStart_;
 739     sliceStarts[3] = openStart_;
 740     sliceStarts[4] = multipleStart_;
 741
 742     forAll(edgeTrees, i)
 743     {
 744         info[i] = edgeTrees[i].findNearest
 745         (
 746             sample,
 747             searchDistSqr[i]
 748         );
 749
 750         // The index returned by the indexedOctree is local to the slice of
 751         // edges it was supplied with, return the index to the value in the
 752         // complete edge list
 753
 754         info[i].setIndex(info[i].index() + sliceStarts[i]);
 755     }
 756 }
 757
 758
 759 const Foam::indexedOctree<Foam::treeDataEdge>&
 760 Foam::extendedFeatureEdgeMesh::edgeTree() const
 761 {
 762     if (edgeTree_.empty())
 763     {
 764         Random rndGen(17301893);
 765
 766         // Slightly extended bb. Slightly off-centred just so on symmetric
 767         // geometry there are less face/edge aligned items.
 768         treeBoundBox bb
 769         (
 770             treeBoundBox(points()).extend(rndGen, 1E-4)
 771         );
 772
 773         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 774         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 775
 776         labelList allEdges(identity(edges().size()));
 777
 778         edgeTree_.reset
 779         (
 780             new indexedOctree<treeDataEdge>
 781             (
 782                 treeDataEdge
 783                 (
 784                     false,          // cachebb
 785                     edges(),        // edges
 786                     points(),       // points
 787                     allEdges        // selected edges
 788                 ),
 789                 bb,     // bb
 790                 8,      // maxLevel
 791                 10,     // leafsize
 792                 3.0     // duplicity
 793             )
 794         );
 795     }
 796
 797     return edgeTree_();
 798 }
 799
 800
 801 const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
 802 Foam::extendedFeatureEdgeMesh::edgeTreesByType() const
 803 {
 804     if (edgeTreesByType_.size() == 0)
 805     {
 806         edgeTreesByType_.setSize(nEdgeTypes);
 807
 808         Random rndGen(872141);
 809
 810         // Slightly extended bb. Slightly off-centred just so on symmetric
 811         // geometry there are less face/edge aligned items.
 812         treeBoundBox bb
 813         (
 814             treeBoundBox(points()).extend(rndGen, 1E-4)
 815         );
 816
 817         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 818         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 819
 820         labelListList sliceEdges(nEdgeTypes);
 821
 822         // External edges
 823         sliceEdges[0] =
 824             identity(internalStart_ - externalStart_) + externalStart_;
 825
 826         // Internal edges
 827         sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
 828
 829         // Flat edges
 830         sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
 831
 832         // Open edges
 833         sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
 834
 835         // Multiple edges
 836         sliceEdges[4] =
 837             identity(edges().size() - multipleStart_) + multipleStart_;
 838
 839         forAll(edgeTreesByType_, i)
 840         {
 841             edgeTreesByType_.set
 842             (
 843                 i,
 844                 new indexedOctree<treeDataEdge>
 845                 (
 846                     treeDataEdge
 847                     (
 848                         false,          // cachebb
 849                         edges(),        // edges
 850                         points(),       // points
 851                         sliceEdges[i]   // selected edges
 852                     ),
 853                     bb,     // bb
 854                     8,      // maxLevel
 855                     10,     // leafsize
 856                     3.0     // duplicity
 857                 )
 858             );
 859         }
 860     }
 861
 862     return edgeTreesByType_;
 863 }
 864
 865
 866 void Foam::extendedFeatureEdgeMesh::writeObj
 867 (
 868     const fileName& prefix
 869 ) const
 870 {
 871     Pout<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
 872         << endl;
 873
 874     label verti = 0;
 875
 876     edgeMesh::write(prefix + "_edgeMesh.obj");
 877
 878     OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
 879     Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
 880
 881     for(label i = 0; i < concaveStart_; i++)
 882     {
 883         meshTools::writeOBJ(convexFtPtStr, points()[i]);
 884     }
 885
 886     OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
 887     Pout<< "Writing concave feature points to "
 888         << concaveFtPtStr.name() << endl;
 889
 890     for(label i = concaveStart_; i < mixedStart_; i++)
 891     {
 892         meshTools::writeOBJ(concaveFtPtStr, points()[i]);
 893     }
 894
 895     OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
 896     Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
 897
 898     for(label i = mixedStart_; i < nonFeatureStart_; i++)
 899     {
 900         meshTools::writeOBJ(mixedFtPtStr, points()[i]);
 901     }
 902
 903     OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
 904     Pout<< "Writing mixed feature point structure to "
 905         << mixedFtPtStructureStr.name() << endl;
 906
 907     verti = 0;
 908     for(label i = mixedStart_; i < nonFeatureStart_; i++)
 909     {
 910         const labelList& ptEds = pointEdges()[i];
 911
 912         forAll(ptEds, j)
 913         {
 914             const edge& e = edges()[ptEds[j]];
 915
 916             meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
 917             meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
 918             mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
 919         }
 920     }
 921
 922     OFstream externalStr(prefix + "_externalEdges.obj");
 923     Pout<< "Writing external edges to " << externalStr.name() << endl;
 924
 925     verti = 0;
 926     for (label i = externalStart_; i < internalStart_; i++)
 927     {
 928         const edge& e = edges()[i];
 929
 930         meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
 931         meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
 932         externalStr << "l " << verti-1 << ' ' << verti << endl;
 933     }
 934
 935     OFstream internalStr(prefix + "_internalEdges.obj");
 936     Pout<< "Writing internal edges to " << internalStr.name() << endl;
 937
 938     verti = 0;
 939     for (label i = internalStart_; i < flatStart_; i++)
 940     {
 941         const edge& e = edges()[i];
 942
 943         meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
 944         meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
 945         internalStr << "l " << verti-1 << ' ' << verti << endl;
 946     }
 947
 948     OFstream flatStr(prefix + "_flatEdges.obj");
 949     Pout<< "Writing flat edges to " << flatStr.name() << endl;
 950
 951     verti = 0;
 952     for (label i = flatStart_; i < openStart_; i++)
 953     {
 954         const edge& e = edges()[i];
 955
 956         meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
 957         meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
 958         flatStr << "l " << verti-1 << ' ' << verti << endl;
 959     }
 960
 961     OFstream openStr(prefix + "_openEdges.obj");
 962     Pout<< "Writing open edges to " << openStr.name() << endl;
 963
 964     verti = 0;
 965     for (label i = openStart_; i < multipleStart_; i++)
 966     {
 967         const edge& e = edges()[i];
 968
 969         meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
 970         meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
 971         openStr << "l " << verti-1 << ' ' << verti << endl;
 972     }
 973
 974     OFstream multipleStr(prefix + "_multipleEdges.obj");
 975     Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
 976
 977     verti = 0;
 978     for (label i = multipleStart_; i < edges().size(); i++)
 979     {
 980         const edge& e = edges()[i];
 981
 982         meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
 983         meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
 984         multipleStr << "l " << verti-1 << ' ' << verti << endl;
 985     }
 986
 987     OFstream regionStr(prefix + "_regionEdges.obj");
 988     Pout<< "Writing region edges to " << regionStr.name() << endl;
 989
 990     verti = 0;
 991     forAll(regionEdges_, i)
 992     {
 993         const edge& e = edges()[regionEdges_[i]];
 994
 995         meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
 996         meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
 997         regionStr << "l " << verti-1 << ' ' << verti << endl;
 998     }
 999 }
1000
1001
1002 bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
1003 {
1004     os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
1005         << "// internalStart, flatStart, openStart, multipleStart, " << nl
1006         << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
1007         << endl;
1008
1009     os  << points() << nl
1010         << edges() << nl
1011         << concaveStart_ << token::SPACE
1012         << mixedStart_ << token::SPACE
1013         << nonFeatureStart_ << token::SPACE
1014         << internalStart_ << token::SPACE
1015         << flatStart_ << token::SPACE
1016         << openStart_ << token::SPACE
1017         << multipleStart_ << nl
1018         << normals_ << nl
1019         << edgeNormals_ << nl
1020         << featurePointNormals_ << nl
1021         << regionEdges_
1022         << endl;
1023
1024     return os.good();
1025 }
1026
1027
1028 // ************************************************************************* //