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) 1991-2008 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Application
  26     surfaceFeatureExtract
  27
  28 Description
  29     Extracts and writes surface features to file
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "triangle.H"
  34 #include "triSurface.H"
  35 #include "argList.H"
  36 #include "surfaceFeatures.H"
  37 #include "treeBoundBox.H"
  38 #include "meshTools.H"
  39 #include "OFstream.H"
  40
  41 using namespace Foam;
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 void dumpBox(const treeBoundBox& bb, const fileName& fName)
  46 {
  47     OFstream str(fName);
  48
  49     Pout<< "Dumping bounding box " << bb << " as lines to obj file "
  50         << str.name() << endl;
  51
  52
  53     pointField boxPoints(bb.points());
  54
  55     forAll(boxPoints, i)
  56     {
  57         meshTools::writeOBJ(str, boxPoints[i]);
  58     }
  59
  60     forAll(treeBoundBox::edges, i)
  61     {
  62         const edge& e = treeBoundBox::edges[i];
  63
  64         str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
  65     }
  66 }
  67
  68
  69 // Deletes all edges inside/outside bounding box from set.
  70 void deleteBox
  71 (
  72     const triSurface& surf,
  73     const treeBoundBox& bb,
  74     const bool removeInside,
  75     List<surfaceFeatures::edgeStatus>& edgeStat
  76 )
  77 {
  78     forAll(edgeStat, edgeI)
  79     {
  80         const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
  81
  82         if
  83         (
  84             (removeInside && bb.contains(eMid))
  85          || (!removeInside && !bb.contains(eMid))
  86         )
  87         {
  88             edgeStat[edgeI] = surfaceFeatures::NONE;
  89         }
  90     }
  91 }
  92
  93
  94 // Main program:
  95
  96 int main(int argc, char *argv[])
  97 {
  98     argList::noParallel();
  99
 100     argList::validArgs.clear();
 101     argList::validArgs.append("surface");
 102     argList::validArgs.append("output set");
 103
 104     argList::validOptions.insert("includedAngle", "included angle [0..180]");
 105     argList::validOptions.insert("set", "input feature set");
 106
 107     argList::validOptions.insert("minLen", "cumulative length of feature");
 108     argList::validOptions.insert("minElem", "number of edges in feature");
 109     argList::validOptions.insert("subsetBox", "((x0 y0 z0)(x1 y1 z1))");
 110     argList::validOptions.insert("deleteBox", "((x0 y0 z0)(x1 y1 z1))");
 111     argList args(argc, argv);
 112
 113     Pout<< "Feature line extraction is only valid on closed manifold surfaces."
 114         << endl;
 115
 116
 117     fileName surfFileName(args.additionalArgs()[0]);
 118     fileName outFileName(args.additionalArgs()[1]);
 119
 120     Pout<< "Surface            : " << surfFileName << nl
 121         << "Output feature set : " << outFileName << nl
 122         << endl;
 123
 124
 125     // Read
 126     // ~~~~
 127
 128     triSurface surf(surfFileName);
 129
 130     Pout<< "Statistics:" << endl;
 131     surf.writeStats(Pout);
 132     Pout<< endl;
 133
 134
 135
 136
 137     // Either construct features from surface&featureangle or read set.
 138     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 139
 140     surfaceFeatures set(surf);
 141
 142     if (args.options().found("set"))
 143     {
 144         fileName setName(args.options()["set"]);
 145
 146         Pout<< "Reading existing feature set from file " << setName << endl;
 147
 148         set = surfaceFeatures(surf, setName);
 149     }
 150     else if (args.options().found("includedAngle"))
 151     {
 152         scalar includedAngle
 153         (
 154             readScalar(IStringStream(args.options()["includedAngle"])())
 155         );
 156
 157         Pout<< "Constructing feature set from included angle " << includedAngle
 158             << endl;
 159
 160         set = surfaceFeatures(surf, includedAngle);
 161
 162         Pout<< endl << "Writing initial features" << endl;
 163         set.write("initial.fSet");
 164         set.writeObj("initial");
 165     }
 166     else
 167     {
 168         FatalErrorIn(args.executable())
 169             << "No initial feature set. Provide either one"
 170             << " of -set (to read existing set)" << nl
 171             << " or -includedAngle (to new set construct from angle)"
 172             << exit(FatalError);
 173     }
 174
 175
 176     Pout<< nl
 177         << "Initial feature set:" << nl
 178         << "    feature points : " << set.featurePoints().size() << nl
 179         << "    feature edges  : " << set.featureEdges().size() << nl
 180         << "    of which" << nl
 181         << "        region edges   : " << set.nRegionEdges() << nl
 182         << "        external edges : " << set.nExternalEdges() << nl
 183         << "        internal edges : " << set.nInternalEdges() << nl
 184         << endl;
 185
 186
 187
 188
 189     // Trim set
 190     // ~~~~~~~~
 191
 192     scalar minLen = -GREAT;
 193     if (args.options().found("minLen"))
 194     {
 195         minLen = readScalar(IStringStream(args.options()["minLen"])());
 196         Pout<< "Removing features of length < " << minLen << endl;
 197     }
 198
 199     label minElem = 0;
 200     if (args.options().found("minElem"))
 201     {
 202         minElem = readLabel(IStringStream(args.options()["minElem"])());
 203         Pout<< "Removing features with number of edges < " << minElem << endl;
 204     }
 205
 206     // Trim away small groups of features
 207     if (minLen > 0 || minLen > 0)
 208     {
 209         set.trimFeatures(minLen, minElem);
 210         Pout<< endl << "Removed small features" << endl;
 211     }
 212
 213
 214
 215     // Subset
 216     // ~~~~~~
 217
 218     // Convert to marked edges, points
 219     List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
 220
 221     if (args.options().found("subsetBox"))
 222     {
 223         treeBoundBox bb(IStringStream(args.options()["subsetBox"])());
 224
 225         Pout<< "Removing all edges outside bb " << bb << endl;
 226         dumpBox(bb, "subsetBox.obj");
 227
 228         deleteBox
 229         (
 230             surf,
 231             bb,
 232             false,
 233             edgeStat
 234         );
 235     }
 236     else if (args.options().found("deleteBox"))
 237     {
 238         treeBoundBox bb(IStringStream(args.options()["deleteBox"])());
 239
 240         Pout<< "Removing all edges inside bb " << bb << endl;
 241         dumpBox(bb, "deleteBox.obj");
 242
 243         deleteBox
 244         (
 245             surf,
 246             bb,
 247             true,
 248             edgeStat
 249         );
 250     }
 251
 252     surfaceFeatures newSet(surf);
 253     newSet.setFromStatus(edgeStat);
 254
 255     Pout<< endl << "Writing trimmed features to " << outFileName << endl;
 256     newSet.write(outFileName);
 257
 258     Pout<< endl << "Writing edge objs." << endl;
 259     newSet.writeObj("final");
 260
 261
 262     Pout<< nl
 263         << "Final feature set:" << nl
 264         << "    feature points : " << newSet.featurePoints().size() << nl
 265         << "    feature edges  : " << newSet.featureEdges().size() << nl
 266         << "    of which" << nl
 267         << "        region edges   : " << newSet.nRegionEdges() << nl
 268         << "        external edges : " << newSet.nExternalEdges() << nl
 269         << "        internal edges : " << newSet.nInternalEdges() << nl
 270         << endl;
 271
 272     Pout<< "End\n" << endl;
 273
 274     return 0;
 275 }
 276
 277
 278 // ************************************************************************* //