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