surfaceHookUp.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-2017 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
27Application
28 surfaceHookUp
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Find close open edges and stitches the surface along them
35
36Usage
37 - surfaceHookUp hookDistance [OPTION]
38
39\*---------------------------------------------------------------------------*/
40
41#include "argList.H"
42#include "Time.H"
43
44#include "triSurfaceMesh.H"
45#include "indexedOctree.H"
46#include "treeBoundBox.H"
47#include "bitSet.H"
48#include "unitConversion.H"
49#include "searchableSurfaces.H"
50#include "IOdictionary.H"
51
52using namespace Foam;
53
54// Split facei along edgeI at position newPointi
55void greenRefine
56(
57 const triSurface& surf,
58 const label facei,
59 const label edgeI,
60 const label newPointi,
62)
63{
64 const labelledTri& f = surf.localFaces()[facei];
65 const edge& e = surf.edges()[edgeI];
66
67 // Find index of edge in face.
68
69 label fp0 = f.find(e[0]);
70 label fp1 = f.fcIndex(fp0);
71 label fp2 = f.fcIndex(fp1);
72
73 if (f[fp1] == e[1])
74 {
75 // Edge oriented like face
76 newFaces.append
77 (
79 (
80 f[fp0],
82 f[fp2],
83 f.region()
84 )
85 );
86 newFaces.append
87 (
89 (
91 f[fp1],
92 f[fp2],
93 f.region()
94 )
95 );
96 }
97 else
98 {
99 newFaces.append
100 (
102 (
103 f[fp2],
104 newPointi,
105 f[fp1],
106 f.region()
107 )
108 );
109 newFaces.append
110 (
112 (
113 newPointi,
114 f[fp0],
115 f[fp1],
116 f.region()
117 )
118 );
119 }
120}
121
122
123//scalar checkEdgeAngle
124//(
125// const triSurface& surf,
126// const label edgeIndex,
127// const label pointIndex,
128// const scalar& angle
129//)
130//{
131// const edge& e = surf.edges()[edgeIndex];
132
133// const vector eVec = e.unitVec(surf.localPoints());
134
135// const labelList& pEdges = surf.pointEdges()[pointIndex];
136//
137// forAll(pEdges, eI)
138// {
139// const edge& nearE = surf.edges()[pEdges[eI]];
140
141// const vector nearEVec = nearE.unitVec(surf.localPoints());
142
143// const scalar dot = eVec & nearEVec;
144// const scalar minCos = degToRad(angle);
145
146// if (mag(dot) > minCos)
147// {
148// return false;
149// }
150// }
151
152// return true;
153//}
154
155
156void createBoundaryEdgeTrees
157(
158 const PtrList<triSurfaceMesh>& surfs,
160 labelListList& treeBoundaryEdges
161)
162{
163 forAll(surfs, surfI)
164 {
165 const triSurface& surf = surfs[surfI];
166
167 // Boundary edges
168 treeBoundaryEdges[surfI] =
170 (
171 surf.nEdges() - surf.nInternalEdges(),
172 surf.nInternalEdges()
173 );
174
175 Random rndGen(17301893);
176
177 // Slightly extended bb. Slightly off-centred just so on symmetric
178 // geometry there are less face/edge aligned items.
179 treeBoundBox bb
180 (
182 );
183
184 bb.min() -= point::uniform(ROOTVSMALL);
185 bb.max() += point::uniform(ROOTVSMALL);
186
187 bEdgeTrees.set
188 (
189 surfI,
191 (
193 (
194 false, // cachebb
195 surf.edges(), // edges
196 surf.localPoints(), // points
197 treeBoundaryEdges[surfI] // selected edges
198 ),
199 bb, // bb
200 8, // maxLevel
201 10, // leafsize
202 3.0 // duplicity
203 )
204 );
205 }
206}
207
208
209class findNearestOpSubset
210{
211 const indexedOctree<treeDataEdge>& tree_;
212
213 DynamicList<label>& shapeMask_;
214
215public:
216
217 findNearestOpSubset
218 (
219 const indexedOctree<treeDataEdge>& tree,
220 DynamicList<label>& shapeMask
221 )
222 :
223 tree_(tree),
224 shapeMask_(shapeMask)
225 {}
226
227 void operator()
228 (
229 const labelUList& indices,
230 const point& sample,
231
232 scalar& nearestDistSqr,
233 label& minIndex,
234 point& nearestPoint
235 ) const
236 {
237 const treeDataEdge& shape = tree_.shapes();
238
239 for (const label index : indices)
240 {
241 const label edgeIndex = shape.edgeLabels()[index];
242
243 if (shapeMask_.found(edgeIndex))
244 {
245 continue;
246 }
247
248 const edge& e = shape.edges()[edgeIndex];
249
250 pointHit nearHit = e.line(shape.points()).nearestDist(sample);
251
252 // Only register hit if closest point is not an edge point
253 if (nearHit.hit())
254 {
255 const scalar distSqr = sqr(nearHit.distance());
256
257 if (distSqr < nearestDistSqr)
258 {
259 nearestDistSqr = distSqr;
260 minIndex = index;
261 nearestPoint = nearHit.rawPoint();
262 }
263 }
264 }
265 }
266};
267
268
269// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
270
271int main(int argc, char *argv[])
272{
273 argList::addNote
274 (
275 "Hook surfaces to other surfaces by moving and retriangulating their"
276 " boundary edges to match other surface boundary edges"
277 );
278 argList::noParallel();
279 argList::addArgument("hookTolerance", "The point merge tolerance");
280 argList::addOption("dict", "file", "Alternative surfaceHookUpDict");
281
282 #include "setRootCase.H"
283 #include "createTime.H"
284
285 const word dictName("surfaceHookUpDict");
287
288 Info<< "Reading " << dictIO.name() << nl << endl;
289
290 const IOdictionary dict(dictIO);
291
292 const scalar dist(args.get<scalar>(1));
293 const scalar matchTolerance(max(1e-6*dist, SMALL));
294 const label maxIters = 100;
295
296 Info<< "Hooking distance = " << dist << endl;
297
299 (
301 (
302 "surfacesToHook",
303 runTime.constant(),
304 "triSurface",
305 runTime
306 ),
307 dict,
308 true // assume single-region names get surface name
309 );
310
311 Info<< nl << "Reading surfaces: " << endl;
312 forAll(surfs, surfI)
313 {
315 Info<< nl << indent << "Surface = " << surfs.names()[surfI] << endl;
316
317 const wordList& regions = surfs[surfI].regions();
318 forAll(regions, surfRegionI)
319 {
321 Info<< indent << "Regions = " << regions[surfRegionI] << endl;
323 }
325 }
326
327 PtrList<indexedOctree<treeDataEdge>> bEdgeTrees(surfs.size());
328 labelListList treeBoundaryEdges(surfs.size());
329
330 List<DynamicList<labelledTri>> newFaces(surfs.size());
331 List<DynamicList<point>> newPoints(surfs.size());
332 List<bitSet> visitedFace(surfs.size());
333
334 PtrList<triSurfaceMesh> newSurfaces(surfs.size());
335 forAll(surfs, surfI)
336 {
337 const triSurfaceMesh& surf =
338 refCast<const triSurfaceMesh>(surfs[surfI]);
339
340 newSurfaces.set
341 (
342 surfI,
344 (
346 (
347 "hookedSurface_" + surfs.names()[surfI],
348 runTime.constant(),
349 "triSurface",
350 runTime
351 ),
352 surf
353 )
354 );
355 }
356
357 label nChanged = 0;
358 label nIters = 1;
359
360 do
361 {
362 Info<< nl << "Iteration = " << nIters++ << endl;
363 nChanged = 0;
364
365 createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);
366
367 forAll(newSurfaces, surfI)
368 {
369 const triSurface& newSurf = newSurfaces[surfI];
370
371 newFaces[surfI] = newSurf.localFaces();
372 newPoints[surfI] = newSurf.localPoints();
373 visitedFace[surfI] = bitSet(newSurf.size(), false);
374 }
375
376 forAll(newSurfaces, surfI)
377 {
378 const triSurface& surf = newSurfaces[surfI];
379
380 List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
381 labelList bPointsHitTree(surf.boundaryPoints().size(), -1);
382
383 const labelListList& pointEdges = surf.pointEdges();
384
385 forAll(bPointsTobEdges, bPointi)
386 {
387 pointIndexHit& nearestHit = bPointsTobEdges[bPointi];
388
389 const label pointi = surf.boundaryPoints()[bPointi];
390 const point& samplePt = surf.localPoints()[pointi];
391
392 const labelList& pEdges = pointEdges[pointi];
393
394 // Add edges connected to the edge to the shapeMask
395 DynamicList<label> shapeMask;
396 shapeMask.append(pEdges);
397
398 forAll(bEdgeTrees, treeI)
399 {
400 const indexedOctree<treeDataEdge>& bEdgeTree =
401 bEdgeTrees[treeI];
402
403 pointIndexHit currentHit =
404 bEdgeTree.findNearest
405 (
406 samplePt,
407 sqr(dist),
408 findNearestOpSubset
409 (
410 bEdgeTree,
411 shapeMask
412 )
413 );
414
415 if
416 (
417 currentHit.hit()
418 &&
419 (
420 !nearestHit.hit()
421 ||
422 (
423 magSqr(currentHit.hitPoint() - samplePt)
424 < magSqr(nearestHit.hitPoint() - samplePt)
425 )
426 )
427 )
428 {
429 nearestHit = currentHit;
430 bPointsHitTree[bPointi] = treeI;
431 }
432 }
433
434 scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);
435
436 if (nearestHit.hit())
437 {
438 // bool rejectEdge =
439 // checkEdgeAngle
440 // (
441 // surf,
442 // nearestHit.index(),
443 // pointi,
444 // 30
445 // );
446
447 if (dist2 > Foam::sqr(dist))
448 {
449 nearestHit.setMiss();
450 }
451 }
452 }
453
454 forAll(bPointsTobEdges, bPointi)
455 {
456 const pointIndexHit& eHit = bPointsTobEdges[bPointi];
457
458 if (eHit.hit())
459 {
460 const label hitSurfI = bPointsHitTree[bPointi];
461 const triSurface& hitSurf = newSurfaces[hitSurfI];
462
463 const label eIndex =
464 treeBoundaryEdges[hitSurfI][eHit.index()];
465 const edge& e = hitSurf.edges()[eIndex];
466
467 const label pointi = surf.boundaryPoints()[bPointi];
468
469 const labelList& eFaces = hitSurf.edgeFaces()[eIndex];
470
471 if (eFaces.size() != 1)
472 {
474 << "Edge is attached to " << eFaces.size()
475 << " faces." << endl;
476
477 continue;
478 }
479
480 const label facei = eFaces[0];
481
482 if (visitedFace[hitSurfI][facei])
483 {
484 continue;
485 }
486
487 DynamicList<labelledTri> newFacesFromSplit(2);
488
489 const point& pt = surf.localPoints()[pointi];
490
491 if
492 (
493 (
494 magSqr(pt - hitSurf.localPoints()[e.start()])
495 < matchTolerance
496 )
497 || (
498 magSqr(pt - hitSurf.localPoints()[e.end()])
499 < matchTolerance
500 )
501 )
502 {
503 continue;
504 }
505
506 nChanged++;
507
508 label newPointi = -1;
509
510 // Keep the points in the same place and move the edge
511 if (hitSurfI == surfI)
512 {
513 newPointi = pointi;
514 }
515 else
516 {
517 newPoints[hitSurfI].append(newPoints[surfI][pointi]);
518 newPointi = newPoints[hitSurfI].size() - 1;
519 }
520
521 // Split the other face.
522 greenRefine
523 (
524 hitSurf,
525 facei,
526 eIndex,
527 newPointi,
528 newFacesFromSplit
529 );
530
531 visitedFace[hitSurfI].set(facei);
532
533 forAll(newFacesFromSplit, newFacei)
534 {
535 const labelledTri& fN = newFacesFromSplit[newFacei];
536
537 if (newFacei == 0)
538 {
539 newFaces[hitSurfI][facei] = fN;
540 }
541 else
542 {
543 newFaces[hitSurfI].append(fN);
544 }
545 }
546 }
547 }
548 }
549
550 Info<< " Number of edges split = " << nChanged << endl;
551
552 forAll(newSurfaces, surfI)
553 {
554 newSurfaces.set
555 (
556 surfI,
558 (
560 (
561 "hookedSurface_" + surfs.names()[surfI],
562 runTime.constant(),
563 "triSurface",
564 runTime
565 ),
567 (
568 newFaces[surfI],
569 newSurfaces[surfI].patches(),
570 pointField(newPoints[surfI])
571 )
572 )
573 );
574 }
575
576 } while (nChanged > 0 && nIters <= maxIters);
577
578 Info<< endl;
579
580 forAll(newSurfaces, surfI)
581 {
582 const triSurfaceMesh& newSurf = newSurfaces[surfI];
583
584 Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
585 << endl;
586
587 newSurf.searchableSurface::write();
588 }
589
590 Info<< "\nEnd\n" << endl;
591
592 return 0;
593}
594
595
596// ************************************************************************* //
Y[inertIndex] max(0.0)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Minimal example by using system/controlDict.functions:
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:54
const point_type & rawPoint() const noexcept
The point, no checks.
Definition: PointHit.H:172
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
const point_type & rawPoint() const noexcept
The point, no checks. Same as point()
label index() const noexcept
Return the hit index.
void setMiss() noexcept
Set the hit status off.
bool hit() const noexcept
Is there a hit?
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
label nEdges() const
Number of edges in patch.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
label nInternalEdges() const
Number of internal edges.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const labelList & boundaryPoints() const
Return list of boundary points, address into LOCAL point list.
const labelListList & edgeFaces() const
Return edge-face addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
Random number generator.
Definition: Random.H:60
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
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
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
const Type & shapes() const
Reference to shape.
A triFace with additional (region) index.
Definition: labelledTri.H:60
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
treeBoundBox extend(Random &rndGen, const scalar s) const
Return slightly wider bounding box.
Holds data for octree to work on an edges subset.
Definition: treeDataEdge.H:57
const edgeList & edges() const
Definition: treeDataEdge.H:173
const labelList & edgeLabels() const
Definition: treeDataEdge.H:183
const pointField & points() const
Definition: treeDataEdge.H:178
IOoject and searching on triSurface.
Triangulated surface description with patch information.
Definition: triSurface.H:79
A class for handling words, derived from Foam::string.
Definition: word.H:68
const polyBoundaryMesh & patches
engineTime & runTime
const word dictName("faMeshDefinition")
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
dimensionedSymmTensor sqr(const dimensionedVector &dv)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
dictionary dict
IOobject dictIO
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.
Random rndGen
Definition: createFields.H:23