FreeFOAM The Cross-Platform CFD Toolkit
surfaceSubset.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  surfaceSubset
26 
27 Description
28  A surface analysis tool which sub-sets the triSurface
29  to choose only a part of interest.
30 
31  Based on subsetMesh.
32 
33 Usage
34 
35  - surfaceSubset [OPTIONS] <surfaceSubsetDic> <Foam surface file> <Foam output surface file>
36 
37  @param <surfaceSubsetDic> \n
38  @todo Detailed description of argument.
39 
40  @param <Foam surface file> \n
41  @todo Detailed description of argument.
42 
43  @param <Foam output surface file> \n
44  @todo Detailed description of argument.
45 
46  @param -case <dir>\n
47  Case directory.
48 
49  @param -help \n
50  Display help message.
51 
52  @param -doc \n
53  Display Doxygen API documentation page for this application.
54 
55  @param -srcDoc \n
56  Display Doxygen source documentation page for this application.
57 
58 \*---------------------------------------------------------------------------*/
59 
60 #include <triSurface/triSurface.H>
61 #include <OpenFOAM/argList.H>
62 #include <OpenFOAM/OFstream.H>
63 #include <OpenFOAM/IFstream.H>
64 #include <OpenFOAM/Switch.H>
65 #include <OpenFOAM/IOdictionary.H>
66 #include <OpenFOAM/boundBox.H>
68 #include <meshTools/octree.H>
70 #include <OpenFOAM/Random.H>
71 
72 using namespace Foam;
73 
74 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
75 // Main program:
76 
77 int main(int argc, char *argv[])
78 {
81  argList::validArgs.append("surfaceSubsetDict");
82  argList::validArgs.append("surface file");
83  argList::validArgs.append("output file");
84  argList args(argc, argv);
85 
86  Info<< "Reading dictionary " << args.additionalArgs()[0] << " ..." << endl;
87  IFstream dictFile(args.additionalArgs()[0]);
88  dictionary meshSubsetDict(dictFile);
89 
90  Info<< "Reading surface " << args.additionalArgs()[1] << " ..." << endl;
91  triSurface surf1(args.additionalArgs()[1]);
92 
93  Info<< "Original:" << endl;
94  surf1.writeStats(Info);
95  Info<< endl;
96 
97 
98  labelList markedPoints
99  (
100  meshSubsetDict.lookup("localPoints")
101  );
102 
103  labelList markedEdges
104  (
105  meshSubsetDict.lookup("edges")
106  );
107 
108  labelList markedFaces
109  (
110  meshSubsetDict.lookup("faces")
111  );
112 
113  pointField markedZone
114  (
115  meshSubsetDict.lookup("zone")
116  );
117 
118  if (markedZone.size() && markedZone.size() != 2)
119  {
121  << "zone specification should be two points, min and max of "
122  << "the boundingbox" << endl
123  << "zone:" << markedZone
124  << exit(FatalError);
125  }
126 
127  Switch addFaceNeighbours
128  (
129  meshSubsetDict.lookup("addFaceNeighbours")
130  );
131 
132  Switch invertSelection
133  (
134  meshSubsetDict.lookup("invertSelection")
135  );
136 
137  // Mark the cells for the subset
138 
139  // Faces to subset
140  boolList facesToSubset(surf1.size(), false);
141 
142 
143  //
144  // pick up faces connected to "localPoints"
145  //
146 
147  if (markedPoints.size())
148  {
149  Info << "Found " << markedPoints.size() << " marked point(s)." << endl;
150 
151  // pick up cells sharing the point
152 
153  forAll (markedPoints, pointI)
154  {
155  if
156  (
157  markedPoints[pointI] < 0
158  || markedPoints[pointI] >= surf1.nPoints()
159  )
160  {
162  << "localPoint label " << markedPoints[pointI]
163  << "out of range."
164  << " The mesh has got "
165  << surf1.nPoints() << " localPoints."
166  << exit(FatalError);
167  }
168 
169  const labelList& curFaces =
170  surf1.pointFaces()[markedPoints[pointI]];
171 
172  forAll (curFaces, i)
173  {
174  facesToSubset[curFaces[i]] = true;
175  }
176  }
177  }
178 
179 
180 
181  //
182  // pick up faces connected to "edges"
183  //
184 
185  if (markedEdges.size())
186  {
187  Info << "Found " << markedEdges.size() << " marked edge(s)." << endl;
188 
189  // pick up cells sharing the edge
190 
191  forAll (markedEdges, edgeI)
192  {
193  if
194  (
195  markedEdges[edgeI] < 0
196  || markedEdges[edgeI] >= surf1.nEdges()
197  )
198  {
200  << "edge label " << markedEdges[edgeI]
201  << "out of range."
202  << " The mesh has got "
203  << surf1.nEdges() << " edges."
204  << exit(FatalError);
205  }
206 
207  const labelList& curFaces = surf1.edgeFaces()[markedEdges[edgeI]];
208 
209  forAll (curFaces, i)
210  {
211  facesToSubset[curFaces[i]] = true;
212  }
213  }
214  }
215 
216 
217  //
218  // pick up faces with centre inside "zone"
219  //
220 
221  if (markedZone.size() == 2)
222  {
223  const point& min = markedZone[0];
224  const point& max = markedZone[1];
225 
226  Info << "Using zone min:" << min << " max:" << max << endl;
227 
228  forAll(surf1, faceI)
229  {
230  const labelledTri& f = surf1[faceI];
231  const point centre = f.centre(surf1.points());
232 
233  if
234  (
235  (centre.x() >= min.x())
236  && (centre.y() >= min.y())
237  && (centre.z() >= min.z())
238  && (centre.x() <= max.x())
239  && (centre.y() <= max.y())
240  && (centre.z() <= max.z())
241  )
242  {
243  facesToSubset[faceI] = true;
244  }
245  }
246  }
247 
248 
249  //
250  // pick up faces on certain side of surface
251  //
252 
253  if (meshSubsetDict.found("surface"))
254  {
255  const dictionary& surfDict = meshSubsetDict.subDict("surface");
256 
257  fileName surfName(surfDict.lookup("name"));
258 
259  Switch outside(surfDict.lookup("outside"));
260 
261  if (outside)
262  {
263  Info<< "Selecting all triangles with centre outside surface "
264  << surfName << endl;
265  }
266  else
267  {
268  Info<< "Selecting all triangles with centre inside surface "
269  << surfName << endl;
270  }
271 
272  // Read surface to select on
273  triSurface selectSurf(surfName);
274 
275  // bb of surface
276  treeBoundBox bb(selectSurf.localPoints());
277 
278  // Radnom number generator
279  Random rndGen(354543);
280 
281  // search engine
283  (
284  treeDataTriSurface(selectSurf),
285  bb.extend(rndGen, 1E-4), // slightly randomize bb
286  8, // maxLevel
287  10, // leafsize
288  3.0 // duplicity
289  );
290 
291  // Check if face (centre) is in outside or inside.
292  forAll(facesToSubset, faceI)
293  {
294  if (!facesToSubset[faceI])
295  {
296  const point fc(surf1[faceI].centre(surf1.points()));
297 
299  selectTree.getVolumeType(fc);
300 
301  if (t == indexedOctree<treeDataTriSurface>::INSIDE && !outside)
302  {
303  facesToSubset[faceI] = true;
304  }
305  else if
306  (
308  && outside
309  )
310  {
311  facesToSubset[faceI] = true;
312  }
313  }
314  }
315  }
316 
317 
318  //
319  // pick up specified "faces"
320  //
321 
322  // Number of additional faces picked up because of addFaceNeighbours
323  label nFaceNeighbours = 0;
324 
325  if (markedFaces.size())
326  {
327  Info << "Found " << markedFaces.size() << " marked face(s)." << endl;
328 
329  // Check and mark faces to pick up
330  forAll (markedFaces, faceI)
331  {
332  if
333  (
334  markedFaces[faceI] < 0
335  || markedFaces[faceI] >= surf1.size()
336  )
337  {
339  << "Face label " << markedFaces[faceI] << "out of range."
340  << " The mesh has got "
341  << surf1.size() << " faces."
342  << exit(FatalError);
343  }
344 
345  // Mark the face
346  facesToSubset[markedFaces[faceI]] = true;
347 
348  // mark its neighbours if requested
349  if (addFaceNeighbours)
350  {
351  const labelList& curFaces =
352  surf1.faceFaces()[markedFaces[faceI]];
353 
354  forAll (curFaces, i)
355  {
356  label faceI = curFaces[i];
357 
358  if (!facesToSubset[faceI])
359  {
360  facesToSubset[faceI] = true;
361  nFaceNeighbours++;
362  }
363  }
364  }
365  }
366  }
367 
368  if (addFaceNeighbours)
369  {
370  Info<< "Added " << nFaceNeighbours
371  << " faces because of addFaceNeighbours" << endl;
372  }
373 
374 
375  if (invertSelection)
376  {
377  Info<< "Inverting selection." << endl;
378  boolList newFacesToSubset(facesToSubset.size());
379 
380  forAll(facesToSubset, i)
381  {
382  if (facesToSubset[i])
383  {
384  newFacesToSubset[i] = false;
385  }
386  else
387  {
388  newFacesToSubset[i] = true;
389  }
390  }
391  facesToSubset.transfer(newFacesToSubset);
392  }
393 
394 
395  // Create subsetted surface
396  labelList pointMap;
397  labelList faceMap;
398  triSurface surf2
399  (
400  surf1.subsetMesh(facesToSubset, pointMap, faceMap)
401  );
402 
403  Info<< "Subset:" << endl;
404  surf2.writeStats(Info);
405  Info << endl;
406 
407  fileName outFileName(args.additionalArgs()[2]);
408 
409  Info << "Writing surface to " << outFileName << endl;
410 
411  surf2.write(outFileName);
412 
413  return 0;
414 }
415 
416 
417 // ************************ vim: set sw=4 sts=4 et: ************************ //