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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 surfaceSubset
29
30Group
31 grpSurfaceUtilities
32
33Description
34 A surface analysis tool that subsets the triSurface to choose a
35 region of interest. Based on subsetMesh.
36
37\*---------------------------------------------------------------------------*/
38
39#include "triSurfaceSearch.H"
40#include "MeshedSurfaces.H"
41#include "argList.H"
42#include "Fstream.H"
43#include "IOdictionary.H"
44#include "boundBox.H"
45#include "indexedOctree.H"
46#include "treeDataTriSurface.H"
47#include "Random.H"
48#include "volumeType.H"
49#include "plane.H"
50
51using namespace Foam;
52
53// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
54
55int main(int argc, char *argv[])
56{
57 argList::addNote
58 (
59 "A surface analysis tool that subsets the surface to choose a"
60 " region of interest."
61 );
62
63 argList::noParallel();
64 argList::addArgument("dict", "The surfaceSubsetDict");
65 argList::addArgument("input", "The input surface file");
66 argList::addArgument("output", "The output surface file");
67 argList args(argc, argv);
68
69 Info<< "Reading dictionary " << args[1] << " ..." << endl;
70 IFstream dictFile(args.get<fileName>(1));
71 dictionary meshSubsetDict(dictFile);
72
73 Info<< "Reading surface " << args[2] << " ..." << endl;
74 meshedSurface surf1(args.get<fileName>(2));
75
76 const auto outFileName(args.get<fileName>(3));
77
78 Info<< "Original:" << endl;
79 surf1.writeStats(Info);
80 Info<< endl;
81
82
83 labelList markedPoints
84 (
85 meshSubsetDict.lookup("localPoints")
86 );
87
88 labelList markedEdges
89 (
90 meshSubsetDict.lookup("edges")
91 );
92
93 labelList markedFaces
94 (
95 meshSubsetDict.lookup("faces")
96 );
97
98 pointField markedZone
99 (
100 meshSubsetDict.lookup("zone")
101 );
102
103
104 boundBox zoneBb;
105
106 if (markedZone.size())
107 {
108 if (markedZone.size() != 2)
109 {
111 << "zone specification should be two points, min and max of "
112 << "the boundingbox" << endl
113 << "zone:" << markedZone
114 << exit(FatalError);
115 }
116
117 zoneBb.min() = markedZone[0];
118 zoneBb.max() = markedZone[1];
119
120 if (!zoneBb.valid())
121 {
123 << "Defined zone is invalid: " << zoneBb << nl;
124 }
125 }
126
127
128 const bool addFaceNeighbours =
129 meshSubsetDict.get<bool>("addFaceNeighbours");
130
131 const bool invertSelection =
132 meshSubsetDict.getOrDefault("invertSelection", false);
133
134 // Mark the cells for the subset
135
136 // Faces to subset
137 bitSet facesToSubset(surf1.size(), false);
138
139
140 //
141 // Faces connected to "localPoints"
142 //
143
144 if (markedPoints.size())
145 {
146 Info<< "Found " << markedPoints.size() << " marked point(s)." << endl;
147
148 for (const label pointi : markedPoints)
149 {
150 if (pointi < 0 || pointi >= surf1.nPoints())
151 {
153 << "localPoint label " << pointi << "out of range."
154 << " Surface has " << surf1.nPoints() << " localPoints."
155 << exit(FatalError);
156 }
157
158 const labelList& curFaces = surf1.pointFaces()[pointi];
159
160 facesToSubset.set(curFaces);
161 }
162 }
163
164
165 //
166 // Faces connected to "edges"
167 //
168
169 if (markedEdges.size())
170 {
171 Info<< "Found " << markedEdges.size() << " marked edge(s)." << endl;
172
173 for (const label edgei : markedEdges)
174 {
175 if (edgei < 0 || edgei >= surf1.nEdges())
176 {
178 << "edge label " << edgei << "out of range."
179 << " Surface has " << surf1.nEdges() << " edges."
180 << exit(FatalError);
181 }
182
183 const labelList& curFaces = surf1.edgeFaces()[edgei];
184
185 facesToSubset.set(curFaces);
186 }
187 }
188
189
190 //
191 // Faces with centre inside "zone"
192 //
193
194 if (zoneBb.valid())
195 {
196 Info<< "Using zone " << zoneBb << endl;
197
198 forAll(surf1, facei)
199 {
200 const point centre = surf1[facei].centre(surf1.points());
201
202 if (zoneBb.contains(centre))
203 {
204 facesToSubset.set(facei);
205 }
206 }
207 }
208
209
210 //
211 // Faces on certain side of surface
212 //
213
214 if (meshSubsetDict.found("surface"))
215 {
216 const dictionary& surfDict = meshSubsetDict.subDict("surface");
217
218 const auto surfName(surfDict.get<fileName>("name"));
219
220 const volumeType::type volType =
221 (
222 surfDict.getOrDefault("outside", false)
223 ? volumeType::OUTSIDE
224 : volumeType::INSIDE
225 );
226
227 Info<< "Selecting faces with centre located "
228 << volumeType::names[volType] << " of surface "
229 << surfName << endl;
230
231 // Read surface to select on
232 triSurface selectSurf(surfName);
233
234 triSurfaceSearch searchSelectSurf
235 (
236 selectSurf,
238 8
239 );
240
241 const indexedOctree<treeDataTriSurface>& selectTree =
242 searchSelectSurf.tree();
243
244 // Check if face (centre) is in outside or inside.
245 forAll(surf1, facei)
246 {
247 if (!facesToSubset[facei])
248 {
249 const point fc(surf1[facei].centre(surf1.points()));
250
251 if (volType == selectTree.getVolumeType(fc))
252 {
253 facesToSubset.set(facei);
254 }
255 }
256 }
257 }
258
259
260 if (meshSubsetDict.found("plane"))
261 {
262 const dictionary& planeDict = meshSubsetDict.subDict("plane");
263
264 const plane pl(planeDict);
265 const scalar distance(planeDict.get<scalar>("distance"));
266 const scalar cosAngle(planeDict.get<scalar>("cosAngle"));
267
268 // Select all triangles that are close to the plane and
269 // whose normal aligns with the plane as well.
270
271 forAll(surf1.faceCentres(), facei)
272 {
273 const point& fc = surf1.faceCentres()[facei];
274 const point& nf = surf1.faceNormals()[facei];
275
276 if (pl.distance(fc) < distance && mag(pl.normal() & nf) > cosAngle)
277 {
278 facesToSubset.set(facei);
279 }
280 }
281 }
282
283
284
285 //
286 // pick up specified "faces"
287 //
288
289 // Number of additional faces picked up because of addFaceNeighbours
290 label nFaceNeighbours = 0;
291
292 if (markedFaces.size())
293 {
294 Info<< "Found " << markedFaces.size() << " marked face(s)." << endl;
295
296 // Check and mark faces to pick up
297 for (const label facei : markedFaces)
298 {
299 if (facei < 0 || facei >= surf1.size())
300 {
302 << "Face label " << facei << "out of range."
303 << " Surface has " << surf1.size() << " faces."
304 << exit(FatalError);
305 }
306
307 // Mark the face
308 facesToSubset.set(facei);
309
310 // Mark its neighbours if requested
311 if (addFaceNeighbours)
312 {
313 const labelList& curFaces = surf1.faceFaces()[facei];
314
315 for (const label neiFacei : curFaces)
316 {
317 if (facesToSubset.set(neiFacei))
318 {
319 ++nFaceNeighbours;
320 }
321 }
322 }
323 }
324 }
325
326 if (addFaceNeighbours)
327 {
328 Info<< "Added " << nFaceNeighbours
329 << " faces because of addFaceNeighbours" << endl;
330 }
331
332
333 if (invertSelection)
334 {
335 Info<< "Inverting selection." << endl;
336
337 facesToSubset.flip();
338 }
339
340
341 // Create subsetted surface
342 meshedSurface surf2(surf1.subsetMesh(facesToSubset));
343
344 Info<< "Subset:" << endl;
345 surf2.writeStats(Info);
346 Info<< endl;
347
348 Info<< "Writing surface to " << outFileName << endl;
349
350 surf2.write(outFileName);
351
352 Info<< "\nEnd\n" << endl;
353
354 return 0;
355}
356
357
358// ************************************************************************* //
Input from file stream, using an ISstream.
Definition: IFstream.H:57
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const Vector< Cmpt > & centre(const Foam::UList< Vector< Cmpt > > &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:114
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
bool valid() const
Bounding box is non-inverted.
Definition: boundBoxI.H:76
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
bool contains(const point &pt) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:271
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
A class for handling file names.
Definition: fileName.H:76
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
Helper class to search on triSurface.
Triangulated surface description with patch information.
Definition: triSurface.H:79
type
Volume classification types.
Definition: volumeType.H:66
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333