cut.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) 2014 OpenFOAM Foundation
9 Copyright (C) 2020 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
27\*---------------------------------------------------------------------------*/
28
29#include "cut.H"
31#include "searchableSurfaces.H"
32#include "triSurfaceMesh.H"
33#include "searchableBox.H"
35#include "surfaceIntersection.H"
36#include "intersectedSurface.H"
37#include "edgeIntersections.H"
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
43namespace searchableSurfaceModifiers
44{
46 addToRunTimeSelectionTable(searchableSurfaceModifier, cut, dictionary);
47}
48}
49
50// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51
53(
54 const faceList& fcs,
55 pointField& pts,
56 triSurface& cutSurf
57) const
58{
59 label nTris = 0;
60 forAll(fcs, i)
61 {
62 nTris += fcs[i].size()-2;
63 }
64
65 DynamicList<labelledTri> tris(nTris);
66
67 forAll(fcs, i)
68 {
69 const face& f = fcs[i];
70 // Triangulate around vertex 0
71 for (label fp = 1; fp < f.size()-1; fp++)
72 {
73 tris.append(labelledTri(f[0], f[fp], f[f.fcIndex(fp)], i));
74 }
75 }
77 forAll(patches, patchi)
78 {
79 patches[patchi] = geometricSurfacePatch
80 (
82 patchi
83 );
84 }
85 cutSurf = triSurface(tris, patches, pts, true);
86}
87
88
90(
91 const searchableSurface& cutter,
92 triSurface& cutSurf
93) const
94{
95 if (isA<searchableBox>(cutter))
96 {
97 const searchableBox& bb = refCast<const searchableBox>(cutter);
98
99 pointField pts(bb.points());
100 triangulate(treeBoundBox::faces, pts, cutSurf);
101
102 return cutSurf;
103 }
104 else if (isA<searchableRotatedBox>(cutter))
105 {
106 const searchableRotatedBox& bb =
107 refCast<const searchableRotatedBox>(cutter);
108
109 pointField pts(bb.points());
110 triangulate(treeBoundBox::faces, pts, cutSurf);
111
112 return cutSurf;
113 }
114 else if (isA<triSurfaceMesh>(cutter))
115 {
116 return const_cast<triSurfaceMesh&>
117 (
118 refCast<const triSurfaceMesh>(cutter)
119 );
120 }
121 else
122 {
124 << "Triangulation only supported for triSurfaceMesh, searchableBox"
125 << ", not for surface " << cutter.name()
126 << " of type " << cutter.type()
127 << exit(FatalError);
128 return const_cast<triSurfaceMesh&>
129 (
130 refCast<const triSurfaceMesh>(cutter)
131 );
132 }
133}
134
135
136// Keep on shuffling surface points until no more degenerate intersections.
137// Moves both surfaces and updates set of edge cuts.
138bool Foam::searchableSurfaceModifiers::cut::intersectSurfaces
139(
140 triSurface& surf1,
141 edgeIntersections& edgeCuts1,
142 triSurface& surf2,
143 edgeIntersections& edgeCuts2
144) const
145{
146 bool hasMoved1 = false;
147 bool hasMoved2 = false;
148
149 for (label iter = 0; iter < 10; iter++)
150 {
151 Info<< "Determining intersections of surf1 edges with surf2"
152 << " faces" << endl;
153
154 // Determine surface1 edge intersections. Allow surface to be moved.
155
156 // Number of iterations needed to resolve degenerates
157 label nIters1 = 0;
158 {
159 triSurfaceSearch querySurf2(surf2);
160
161 scalarField surf1PointTol
162 (
164 );
165
166 // Determine raw intersections
167 edgeCuts1 = edgeIntersections
168 (
169 surf1,
170 querySurf2,
171 surf1PointTol
172 );
173
174 // Shuffle a bit to resolve degenerate edge-face hits
175 {
176 pointField points1(surf1.points());
177
178 nIters1 =
179 edgeCuts1.removeDegenerates
180 (
181 5, // max iterations
182 surf1,
183 querySurf2,
184 surf1PointTol,
185 points1 // work array
186 );
187
188 if (nIters1 != 0)
189 {
190 // Update geometric quantities
191 surf1.movePoints(points1);
192 hasMoved1 = true;
193 }
194 }
195 }
196
197 Info<< "Determining intersections of surf2 edges with surf1"
198 << " faces" << endl;
199
200 label nIters2 = 0;
201 {
202 triSurfaceSearch querySurf1(surf1);
203
204 scalarField surf2PointTol
205 (
207 );
208
209 // Determine raw intersections
210 edgeCuts2 = edgeIntersections
211 (
212 surf2,
213 querySurf1,
214 surf2PointTol
215 );
216
217 // Shuffle a bit to resolve degenerate edge-face hits
218 {
219 pointField points2(surf2.points());
220
221 nIters2 =
222 edgeCuts2.removeDegenerates
223 (
224 5, // max iterations
225 surf2,
226 querySurf1,
227 surf2PointTol,
228 points2 // work array
229 );
230
231 if (nIters2 != 0)
232 {
233 // Update geometric quantities
234 surf2.movePoints(points2);
235 hasMoved2 = true;
236 }
237 }
238 }
239
240
241 if (nIters1 == 0 && nIters2 == 0)
242 {
243 //Info<< "** Resolved all intersections to be proper"
244 // << "edge-face pierce" << endl;
245 break;
246 }
247 }
248
249 //if (hasMoved1)
250 //{
251 // fileName newFile("surf1.ftr");
252 // Info<< "Surface 1 has been moved. Writing to " << newFile
253 // << endl;
254 // surf1.write(newFile);
255 //}
256 //
257 //if (hasMoved2)
258 //{
259 // fileName newFile("surf2.ftr");
260 // Info<< "Surface 2 has been moved. Writing to " << newFile
261 // << endl;
262 // surf2.write(newFile);
263 //}
264
265 return hasMoved1 || hasMoved2;
266}
267
268
269// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
270
272(
273 const searchableSurfaces& geometry,
274 const dictionary& dict
275)
276:
277 searchableSurfaceModifier(geometry, dict),
278 cutterNames_(dict_.get<wordRes>("cutters"))
279{}
280
281
282// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
283
285(
286 const labelList& regions,
287 searchableSurface& geom
288) const
289{
290 triSurface& surf = refCast<triSurfaceMesh>(geom);
291
292 bool changed = false;
293
294 // Find the surfaces to cut with
295 for (const wordRe& cutterName : cutterNames_)
296 {
297 labelList geomIDs = findStrings(cutterName, geometry_.names());
298
299 for (const label geomI : geomIDs)
300 {
301 const searchableSurface& cutter = geometry_[geomI];
302
303 // Triangulate
304 triSurface work;
305 triSurface& cutSurf = triangulate(cutter, work);
306
307 // Determine intersection (with perturbation)
308 edgeIntersections edge1Cuts;
309 edgeIntersections edge2Cuts;
310 intersectSurfaces
311 (
312 surf,
313 edge1Cuts,
314 cutSurf,
315 edge2Cuts
316 );
317
318
319 // Determine intersection edges
320 surfaceIntersection inter(surf, edge1Cuts, cutSurf, edge2Cuts);
321
322
323 // Use intersection edges to cut up faces. (does all the hard work)
324 intersectedSurface surf3(surf, true, inter);
325
326
327 // Mark triangles based on whether they are inside or outside
328 List<volumeType> volTypes;
329 cutter.getVolumeType(surf3.faceCentres(), volTypes);
330
331 label nInside = 0;
332 forAll(volTypes, i)
333 {
334 if (volTypes[i] == volumeType::INSIDE)
335 {
336 ++nInside;
337 }
338 }
339
340 // Add a patch for inside the box
341 if (nInside && surf3.patches().size() > 0)
342 {
343 geometricSurfacePatchList newPatches(surf3.patches());
344 label sz = newPatches.size();
345 newPatches.setSize(sz+1);
346 newPatches[sz] = geometricSurfacePatch
347 (
348 newPatches[sz-1].name() + "_inside",
349 newPatches[sz-1].index(),
350 newPatches[sz-1].geometricType()
351 );
352
353 Info<< "Moving " << nInside << " out of " << surf3.size()
354 << " triangles to region "
355 << newPatches[sz].name() << endl;
356
357
358 List<labelledTri> newTris(surf3);
359 forAll(volTypes, i)
360 {
361 if (volTypes[i] == volumeType::INSIDE)
362 {
363 newTris[i].region() = sz;
364 }
365 }
366 pointField newPoints(surf3.points());
367 surf = triSurface(newTris, newPatches, newPoints, true);
368
369 changed = true;
370 }
371 }
372 }
373
374 return changed;
375}
376
377
378// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
virtual label triangulate()
Triangulate in-place, returning the number of triangles added.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
static word defaultName
The default cloud name: defaultCloud.
Definition: cloud.H:90
static scalarField minEdgeLength(const triSurface &surf)
Calculate min edge length for every surface point.
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
const searchableSurfaces & geometry_
virtual bool modify(const labelList &regions, searchableSurface &) const
Apply this selector.
const wordList & names() const
Surface names, not region names.
static const faceList faces
Face to point addressing.
Definition: treeBoundBox.H:151
Triangulated surface description with patch information.
Definition: triSurface.H:79
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
List< geometricSurfacePatch > geometricSurfacePatchList
A List of geometricSurfacePatch.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
labelList findStrings(const regExp &matcher, const UList< StringType > &input, const bool invert=false)
Return list indices for strings matching the regular expression.
Definition: stringListOps.H:87
labelList f(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333