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