applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.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 Application
  25     surfaceFeatureExtract
  26
  27 Description
  28     Extracts and writes surface features to file
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32
  33 #include "triangle.H"
  34 #include "triSurface.H"
  35 #include "argList.H"
  36 #include "Time.H"
  37 #include "surfaceFeatures.H"
  38 #include "featureEdgeMesh.H"
  39 #include "extendedFeatureEdgeMesh.H"
  40 #include "treeBoundBox.H"
  41 #include "meshTools.H"
  42 #include "OFstream.H"
  43 #include "unitConversion.H"
  44
  45 using namespace Foam;
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 void dumpBox(const treeBoundBox& bb, const fileName& fName)
  50 {
  51     OFstream str(fName);
  52
  53     Info<< "Dumping bounding box " << bb << " as lines to obj file "
  54         << str.name() << endl;
  55
  56
  57     pointField boxPoints(bb.points());
  58
  59     forAll(boxPoints, i)
  60     {
  61         meshTools::writeOBJ(str, boxPoints[i]);
  62     }
  63
  64     forAll(treeBoundBox::edges, i)
  65     {
  66         const edge& e = treeBoundBox::edges[i];
  67
  68         str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
  69     }
  70 }
  71
  72
  73 // Deletes all edges inside/outside bounding box from set.
  74 void deleteBox
  75 (
  76     const triSurface& surf,
  77     const treeBoundBox& bb,
  78     const bool removeInside,
  79     List<surfaceFeatures::edgeStatus>& edgeStat
  80 )
  81 {
  82     forAll(edgeStat, edgeI)
  83     {
  84         const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
  85
  86         if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
  87         {
  88             edgeStat[edgeI] = surfaceFeatures::NONE;
  89         }
  90     }
  91 }
  92
  93
  94 // Unmark non-manifold edges if individual triangles are not features
  95 void unmarkBaffles
  96 (
  97     const triSurface& surf,
  98     const scalar includedAngle,
  99     List<surfaceFeatures::edgeStatus>& edgeStat
 100 )
 101 {
 102     scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
 103
 104     const labelListList& edgeFaces = surf.edgeFaces();
 105
 106     forAll(edgeFaces, edgeI)
 107     {
 108         const labelList& eFaces = edgeFaces[edgeI];
 109
 110         if (eFaces.size() > 2)
 111         {
 112             label i0 = eFaces[0];
 113             //const labelledTri& f0 = surf[i0];
 114             const vector& n0 = surf.faceNormals()[i0];
 115
 116             //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
 117
 118             bool same = true;
 119
 120             for (label i = 1; i < eFaces.size(); i++)
 121             {
 122                 //const labelledTri& f = surf[i];
 123                 const vector& n = surf.faceNormals()[eFaces[i]];
 124
 125                 //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
 126
 127                 if (mag(n&n0) < minCos)
 128                 {
 129                     same = false;
 130                     break;
 131                 }
 132             }
 133
 134             if (same)
 135             {
 136                 edgeStat[edgeI] = surfaceFeatures::NONE;
 137             }
 138         }
 139     }
 140 }
 141
 142
 143 // Main program:
 144
 145 int main(int argc, char *argv[])
 146 {
 147     argList::addNote
 148     (
 149         "extract and write surface features to file"
 150     );
 151     argList::noParallel();
 152     argList::validArgs.append("surface");
 153     argList::validArgs.append("output set");
 154
 155     argList::addOption
 156     (
 157         "includedAngle",
 158         "degrees",
 159         "construct feature set from included angle [0..180]"
 160     );
 161     argList::addOption
 162     (
 163         "set",
 164         "name",
 165         "use existing feature set from file"
 166     );
 167     argList::addOption
 168     (
 169         "minLen",
 170         "scalar",
 171         "remove features shorter than the specified cumulative length"
 172     );
 173     argList::addOption
 174     (
 175         "minElem",
 176         "int",
 177         "remove features with fewer than the specified number of edges"
 178     );
 179     argList::addOption
 180     (
 181         "subsetBox",
 182         "((x0 y0 z0)(x1 y1 z1))",
 183         "remove edges outside specified bounding box"
 184     );
 185     argList::addOption
 186     (
 187         "deleteBox",
 188         "((x0 y0 z0)(x1 y1 z1))",
 189         "remove edges within specified bounding box"
 190     );
 191     argList::addBoolOption
 192     (
 193         "writeObj",
 194         "write extendedFeatureEdgeMesh obj files"
 195     );
 196
 197 #   include "setRootCase.H"
 198 #   include "createTime.H"
 199
 200     Info<< "Feature line extraction is only valid on closed manifold surfaces."
 201         << endl;
 202
 203     bool writeObj = args.optionFound("writeObj");
 204
 205     const fileName surfFileName = args[1];
 206     const fileName outFileName  = args[2];
 207
 208     Info<< "Surface            : " << surfFileName << nl
 209         << "Output feature set : " << outFileName << nl
 210         << endl;
 211
 212     fileName sFeatFileName = surfFileName.lessExt().name();
 213
 214
 215     // Read
 216     // ~~~~
 217
 218     triSurface surf(surfFileName);
 219
 220     Info<< "Statistics:" << endl;
 221     surf.writeStats(Info);
 222     Info<< endl;
 223
 224
 225     // Either construct features from surface&featureangle or read set.
 226     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 227
 228     surfaceFeatures set(surf);
 229
 230     if (args.optionFound("set"))
 231     {
 232         const fileName setName = args["set"];
 233
 234         Info<< "Reading existing feature set from file " << setName << endl;
 235
 236         set = surfaceFeatures(surf, setName);
 237     }
 238     else if (args.optionFound("includedAngle"))
 239     {
 240         const scalar includedAngle = args.optionRead<scalar>("includedAngle");
 241
 242         Info<< "Constructing feature set from included angle " << includedAngle
 243             << endl;
 244
 245         set = surfaceFeatures(surf, includedAngle);
 246
 247         // Info<< nl << "Writing initial features" << endl;
 248         // set.write("initial.fSet");
 249         // set.writeObj("initial");
 250     }
 251     else
 252     {
 253         FatalErrorIn(args.executable())
 254             << "No initial feature set. Provide either one"
 255             << " of -set (to read existing set)" << nl
 256             << " or -includedAngle (to new set construct from angle)"
 257             << exit(FatalError);
 258     }
 259
 260     Info<< nl
 261         << "Initial feature set:" << nl
 262         << "    feature points : " << set.featurePoints().size() << nl
 263         << "    feature edges  : " << set.featureEdges().size() << nl
 264         << "    of which" << nl
 265         << "        region edges   : " << set.nRegionEdges() << nl
 266         << "        external edges : " << set.nExternalEdges() << nl
 267         << "        internal edges : " << set.nInternalEdges() << nl
 268         << endl;
 269
 270
 271     // Trim set
 272     // ~~~~~~~~
 273
 274     scalar minLen = -GREAT;
 275     if (args.optionReadIfPresent("minLen", minLen))
 276     {
 277         Info<< "Removing features of length < " << minLen << endl;
 278     }
 279
 280     label minElem = 0;
 281     if (args.optionReadIfPresent("minElem", minElem))
 282     {
 283         Info<< "Removing features with number of edges < " << minElem << endl;
 284     }
 285
 286     // Trim away small groups of features
 287     if (minElem > 0 || minLen > 0)
 288     {
 289         set.trimFeatures(minLen, minElem);
 290         Info<< endl << "Removed small features" << endl;
 291     }
 292
 293
 294     // Subset
 295     // ~~~~~~
 296
 297     // Convert to marked edges, points
 298     List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
 299
 300     if (args.optionFound("subsetBox"))
 301     {
 302         treeBoundBox bb
 303         (
 304             args.optionLookup("subsetBox")()
 305         );
 306
 307         Info<< "Removing all edges outside bb " << bb << endl;
 308         dumpBox(bb, "subsetBox.obj");
 309
 310         deleteBox
 311         (
 312             surf,
 313             bb,
 314             false,
 315             edgeStat
 316         );
 317     }
 318     else if (args.optionFound("deleteBox"))
 319     {
 320         treeBoundBox bb
 321         (
 322             args.optionLookup("deleteBox")()
 323         );
 324
 325         Info<< "Removing all edges inside bb " << bb << endl;
 326         dumpBox(bb, "deleteBox.obj");
 327
 328         deleteBox
 329         (
 330             surf,
 331             bb,
 332             true,
 333             edgeStat
 334         );
 335     }
 336
 337     surfaceFeatures newSet(surf);
 338     newSet.setFromStatus(edgeStat);
 339
 340     Info<< endl << "Writing trimmed features to " << outFileName << endl;
 341     newSet.write(outFileName);
 342
 343     // Info<< endl << "Writing edge objs." << endl;
 344     // newSet.writeObj("final");
 345
 346     Info<< nl
 347         << "Final feature set:" << nl
 348         << "    feature points : " << newSet.featurePoints().size() << nl
 349         << "    feature edges  : " << newSet.featureEdges().size() << nl
 350         << "    of which" << nl
 351         << "        region edges   : " << newSet.nRegionEdges() << nl
 352         << "        external edges : " << newSet.nExternalEdges() << nl
 353         << "        internal edges : " << newSet.nInternalEdges() << nl
 354         << endl;
 355
 356     // Extracting and writing a extendedFeatureEdgeMesh
 357
 358     extendedFeatureEdgeMesh feMesh
 359     (
 360         newSet,
 361         runTime,
 362         sFeatFileName + ".extendedFeatureEdgeMesh"
 363     );
 364
 365     Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
 366         << endl;
 367
 368     if (writeObj)
 369     {
 370         feMesh.writeObj(surfFileName.lessExt().name());
 371     }
 372
 373     feMesh.write();
 374
 375
 376     // Write a featureEdgeMesh for backwards compatibility
 377     {
 378         featureEdgeMesh bfeMesh
 379         (
 380             IOobject
 381             (
 382                 surfFileName.lessExt().name() + ".eMesh",   // name
 383                 runTime.constant(), // instance
 384                 "triSurface",
 385                 runTime,            // registry
 386                 IOobject::NO_READ,
 387                 IOobject::AUTO_WRITE,
 388                 false
 389             ),
 390             feMesh.points(),
 391             feMesh.edges()
 392         );
 393
 394         Info<< nl << "Writing featureEdgeMesh to "
 395             << bfeMesh.objectPath() << endl;
 396
 397         bfeMesh.regIOobject::write();
 398     }
 399
 400     Info<< "End\n" << endl;
 401
 402     return 0;
 403 }
 404
 405
 406 // ************************************************************************* //