foamyHexMeshSurfaceSimplify.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) 2012-2016 OpenFOAM Foundation
9 Copyright (C) 2019-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 foamyHexMeshSurfaceSimplify
29
30Description
31 Simplifies surfaces by resampling.
32
33 Uses Thomas Lewiner's topology preserving MarchingCubes.
34 (http://zeus.mat.puc-rio.br/tomlew/tomlew_uk.php)
35
36\*---------------------------------------------------------------------------*/
37
38#include "argList.H"
39#include "Time.H"
40#include "searchableSurfaces.H"
42#include "triSurfaceMesh.H"
43
44#include "opt_octree.h"
45#include "cube.h"
46
47using namespace Foam;
48
49// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50
52{
53 const vector scale_;
54
55 const vector offset_;
56
57public:
58
59 //- Construct from components
61 (
62 const vector scale,
63 const vector offset
64 )
65 :
66 scale_(scale),
67 offset_(offset)
68 {}
69
70 inline Point toLocal(const Foam::point& pt) const
71 {
72 Foam::point p = cmptMultiply(scale_, (pt + offset_));
73 return Point(p.x(), p.y(), p.z());
74 }
75
76 inline Foam::point toGlobal(const Point& pt) const
77 {
78 point p(pt.x(), pt.y(), pt.z());
79 return cmptDivide(p, scale_) - offset_;
80 }
81};
82
83
84
85
86// For use in Fast-Dual Octree from Thomas Lewiner
87class distanceCalc
88:
89 public ::data_access
90{
91
92 const Level min_level_;
93
94 const conformationSurfaces& geometryToConformTo_;
95
96 const pointConversion& converter_;
97
98
99 // Private Member Functions
100
101 scalar signedDistance(const Foam::point& pt) const
102 {
103 const searchableSurfaces& geometry = geometryToConformTo_.geometry();
104 const labelList& surfaces = geometryToConformTo_.surfaces();
105
106 static labelList nearestSurfaces;
107 static scalarField distance;
108
109 static pointField samples(1);
110 samples[0] = pt;
111
112 searchableSurfacesQueries::signedDistance
113 (
114 geometry,
115 surfaces,
116 samples,
117 scalarField(1, GREAT),
118 volumeType::OUTSIDE,
119 nearestSurfaces,
120 distance
121 );
122
123 return distance[0];
124 }
125
126
127public:
128
129 // Constructors
130
131 //- Construct from components
132 distanceCalc
133 (
134 Level max_level_,
135 real iso_val_,
136 Level min_level,
137 const conformationSurfaces& geometryToConformTo,
138 const pointConversion& converter
139 )
140 :
141 data_access(max_level_,iso_val_),
142 min_level_(min_level),
143 geometryToConformTo_(geometryToConformTo),
144 converter_(converter)
145 {}
146
147
148 //- Destructor
149 virtual ~distanceCalc() = default;
150
151
152 // Member Functions
153
154 //- Test function
155 virtual bool need_refine( const Cube &c )
156 {
157 int l = c.lv() ;
158
159 if ( l >= _max_level ) return false;
160 if ( l < min_level_ ) return true;
161
162 treeBoundBox bb
163 (
164 converter_.toGlobal
165 (
166 Point
167 (
168 c.xmin(),
169 c.ymin(),
170 c.zmin()
171 )
172 ),
173 converter_.toGlobal
174 (
175 Point
176 (
177 c.xmax(),
178 c.ymax(),
179 c.zmax()
180 )
181 )
182 );
183
184 const searchableSurfaces& geometry =
185 geometryToConformTo_.geometry();
186 const labelList& surfaces =
187 geometryToConformTo_.surfaces();
188
189
190 //- Uniform refinement around surface
191 {
192 forAll(surfaces, i)
193 {
194 if (geometry[surfaces[i]].overlaps(bb))
195 {
196 return true;
197 }
198 }
199 return false;
200 }
201
202
204 //scalar ccDist = signedDistance(bb.centre());
205 //scalar ccVal = ccDist - _iso_val;
206 //if (mag(ccVal) < SMALL)
207 //{
208 // return true;
209 //}
210 //const pointField points(bb.points());
211 //forAll(points, pointi)
212 //{
213 // scalar pointVal = signedDistance(points[pointi]) - _iso_val;
214 // if (ccVal*pointVal < 0)
215 // {
216 // return true;
217 // }
218 //}
219 //return false;
220
221
225 //const pointField points(bb.points());
226 //const edgeList& edges = treeBoundBox::edges;
227 //pointField start(edges.size());
228 //pointField end(edges.size());
229 //forAll(edges, i)
230 //{
231 // start[i] = points[edges[i][0]];
232 // end[i] = points[edges[i][1]];
233 //}
234 //Foam::List<Foam::List<pointIndexHit>> hitInfo;
235 //labelListList hitSurfaces;
236 //searchableSurfacesQueries::findAllIntersections
237 //(
238 // geometry,
239 // surfaces,
240 // start,
241 // end,
242 // hitSurfaces,
243 // hitInfo
244 //);
245 //
247 //label nInt = 0;
248 //forAll(hitSurfaces, edgeI)
249 //{
250 // nInt += hitSurfaces[edgeI].size();
251 //}
252 //
253 //if (nInt == 0)
254 //{
255 // // No surface intersected. See if there is one inside
256 // forAll(surfaces, i)
257 // {
258 // if (geometry[surfaces[i]].overlaps(bb))
259 // {
260 // return true;
261 // }
262 // }
263 // return false;
264 //}
265 //
267 //label baseSurfI = -1;
268 //forAll(hitSurfaces, edgeI)
269 //{
270 // const labelList& hSurfs = hitSurfaces[edgeI];
271 // forAll(hSurfs, i)
272 // {
273 // if (baseSurfI == -1)
274 // {
275 // baseSurfI = hSurfs[i];
276 // }
277 // else if (baseSurfI != hSurfs[i])
278 // {
279 // // Multiple surfaces
280 // return true;
281 // }
282 // }
283 //}
284 //
286 //DynamicList<pointIndexHit> baseInfo(nInt);
287 //forAll(hitInfo, edgeI)
288 //{
289 // const Foam::List<pointIndexHit>& hits = hitInfo[edgeI];
290 // forAll(hits, i)
291 // {
292 // (void)hits[i].hitPoint();
293 // baseInfo.append(hits[i]);
294 // }
295 //}
296 //vectorField normals;
297 //geometry[surfaces[baseSurfI]].getNormal(baseInfo, normals);
298 //for (label i = 1; i < normals.size(); ++i)
299 //{
300 // if ((normals[0] & normals[i]) < 0.9)
301 // {
302 // return true;
303 // }
304 //}
305 //labelList regions;
306 //geometry[surfaces[baseSurfI]].getRegion(baseInfo, regions);
307 //for (label i = 1; i < regions.size(); ++i)
308 //{
309 // if (regions[0] != regions[i])
310 // {
311 // return true;
312 // }
313 //}
314 //return false;
315
316
317
318 //samples[0] = point(c.xmin(), c.ymin(), c.zmin());
319 //samples[1] = point(c.xmax(), c.ymin(), c.zmin());
320 //samples[2] = point(c.xmax(), c.ymax(), c.zmin());
321 //samples[3] = point(c.xmin(), c.ymax(), c.zmin());
322 //
323 //samples[4] = point(c.xmin(), c.ymin(), c.zmax());
324 //samples[5] = point(c.xmax(), c.ymin(), c.zmax());
325 //samples[6] = point(c.xmax(), c.ymax(), c.zmax());
326 //samples[7] = point(c.xmin(), c.ymax(), c.zmax());
327
328 //scalarField nearestDistSqr(8, GREAT);
329 //
330 //Foam::List<pointIndexHit> nearestInfo;
331 //surf_.findNearest(samples, nearestDistSqr, nearestInfo);
332 //vectorField normals;
333 //surf_.getNormal(nearestInfo, normals);
334 //
335 //for (label i = 1; i < normals.size(); ++i)
336 //{
337 // if ((normals[0] & normals[i]) < 0.5)
338 // {
339 // return true;
340 // }
341 //}
342 //return false;
343
345 //const labelList elems(surf_.tree().findBox(bb));
346 //
347 //if (elems.size() > 1)
348 //{
349 // return true;
350 //}
351 //else
352 //{
353 // return false;
354 //}
355 }
356
357 //- Data function
358 virtual real value_at( const Cube &c )
359 {
360 return signedDistance(converter_.toGlobal(c)) - _iso_val;
361 }
362};
363
364
365// Main program:
366
367int main(int argc, char *argv[])
368{
369 argList::addNote
370 (
371 "Re-sample surfaces used in foamyHexMesh operation"
372 );
373
374 argList::addArgument("output", "The output triSurface/ file");
375
376 argList::noFunctionObjects(); // Never use function objects
377
378 #include "setRootCase.H"
379 #include "createTime.H"
380
381 const auto exportName = args.get<fileName>(1);
382
383 Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
384 << " writing a re-sampled surface to " << exportName
385 << nl << endl;
386
388
389 IOdictionary foamyHexMeshDict
390 (
392 (
393 "foamyHexMeshDict",
394 runTime.system(),
395 runTime,
396 IOobject::MUST_READ_IF_MODIFIED,
397 IOobject::NO_WRITE
398 )
399 );
400
401 // Define/load all geometry
402 searchableSurfaces allGeometry
403 (
405 (
406 "cvSearchableSurfaces",
407 runTime.constant(),
408 "triSurface",
409 runTime,
410 IOobject::MUST_READ,
411 IOobject::NO_WRITE
412 ),
413 foamyHexMeshDict.subDict("geometry"),
414 foamyHexMeshDict.getOrDefault("singleRegionName", true)
415 );
416
417 Info<< "Geometry read in = "
418 << timer.cpuTimeIncrement() << " s." << nl << endl;
419
420
421 Random rndGen(64293*Pstream::myProcNo());
422
423 conformationSurfaces geometryToConformTo
424 (
425 runTime,
426 rndGen,
427 allGeometry,
428 foamyHexMeshDict.subDict("surfaceConformation")
429 );
430
431 Info<< "Set up geometry in = "
432 << timer.cpuTimeIncrement() << " s." << nl << endl;
433
434
435 const searchableSurfaces& geometry = geometryToConformTo.geometry();
436 const labelList& surfaces = geometryToConformTo.surfaces();
437
438
439 const label minLevel = 2;
440
441 // The max cube size follows from the minLevel and the default cube size
442 // (1)
443 const scalar maxSize = 1.0 / (1 << minLevel);
444 const scalar halfMaxSize = 0.5*maxSize;
445
446
447 // Scale the geometry to fit within
448 // halfMaxSize .. 1-halfMaxSize
449
450 scalar wantedRange = 1.0-maxSize;
451
452 const treeBoundBox bb = geometryToConformTo.globalBounds();
453
454 const vector scale = cmptDivide
455 (
456 point(wantedRange, wantedRange, wantedRange),
457 bb.span()
458 );
459 const vector offset =
461 (
462 point(halfMaxSize, halfMaxSize, halfMaxSize),
463 scale
464 )
465 -bb.min();
466
467
468 const pointConversion converter(scale, offset);
469
470
471 // Marching cubes
472
473 OptOctree octree;
474
475 distanceCalc ref
476 (
477 8, //maxLevel
478 0.0, //distance
479 minLevel, //minLevel
480 geometryToConformTo,
481 converter
482 );
483
484 octree.refine(&ref);
485 octree.set_impl(&ref);
486
487 Info<< "Calculated octree in = "
488 << timer.cpuTimeIncrement() << " s." << nl << endl;
489
490 MarchingCubes& mc = octree.mc();
491
492 mc.clean_all() ;
493 octree.build_isosurface(&ref) ;
494
495 Info<< "Constructed iso surface of distance in = "
496 << timer.cpuTimeIncrement() << " s." << nl << endl;
497
498 // Write output file
499 if (mc.ntrigs() > 0)
500 {
501 Triangle* triangles = mc.triangles();
502 label nTris = mc.ntrigs();
503 Foam::DynamicList<labelledTri> tris(mc.ntrigs());
504 for (label triI = 0; triI < nTris; ++triI)
505 {
506 const Triangle& t = triangles[triI];
507 if (t.v1 != t.v2 && t.v1 != t.v3 && t.v2 != t.v3)
508 {
509 tris.append
510 (
512 (
513 triangles[triI].v1,
514 triangles[triI].v2,
515 triangles[triI].v3,
516 0 // region
517 )
518 );
519 }
520 }
521
522
523 Point* vertices = mc.vertices();
524 pointField points(mc.nverts());
525 forAll(points, pointi)
526 {
527 const Point& v = vertices[pointi];
528 points[pointi] = converter.toGlobal(v);
529 }
530
531
532 // Find correspondence to original surfaces
533 labelList regionOffsets(surfaces.size());
534 label nRegions = 0;
535 forAll(surfaces, i)
536 {
537 const wordList& regions = geometry[surfaces[i]].regions();
538 regionOffsets[i] = nRegions;
539 nRegions += regions.size();
540 }
541
542
544 nRegions = 0;
545 forAll(surfaces, i)
546 {
547 const wordList& regions = geometry[surfaces[i]].regions();
548
549 forAll(regions, regionI)
550 {
552 (
553 geometry[surfaces[i]].name() + "_" + regions[regionI],
554 nRegions,
555 "patch"
556 );
557 nRegions++;
558 }
559 }
560
561 triSurface s(tris, patches, points, true);
562 tris.clearStorage();
563
564 Info<< "Extracted triSurface in = "
565 << timer.cpuTimeIncrement() << " s." << nl << endl;
566
567
568 // Find out region on local surface of nearest point
569 {
571 labelList hitSurfaces;
572 geometryToConformTo.findSurfaceNearest
573 (
574 s.faceCentres(),
575 scalarField(s.size(), sqr(GREAT)),
576 hitInfo,
577 hitSurfaces
578 );
579
580 // Get region
581 DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
582 DynamicList<label> surfIndices(hitSurfaces.size());
583
584 forAll(surfaces, surfI)
585 {
586 // Extract info on this surface
587 surfInfo.clear();
588 surfIndices.clear();
589 forAll(hitSurfaces, triI)
590 {
591 if (hitSurfaces[triI] == surfI)
592 {
593 surfInfo.append(hitInfo[triI]);
594 surfIndices.append(triI);
595 }
596 }
597
598 // Calculate sideness of these surface points
599 labelList region;
600 geometry[surfaces[surfI]].getRegion(surfInfo, region);
601
602 forAll(region, i)
603 {
604 label triI = surfIndices[i];
605 s[triI].region() = regionOffsets[surfI]+region[i];
606 }
607 }
608 }
609
610 Info<< "Re-patched surface in = "
611 << timer.cpuTimeIncrement() << " s." << nl << endl;
612
613 triSurfaceMesh smesh
614 (
616 (
617 exportName,
618 runTime.constant(), // instance
619 "triSurface",
620 runTime, // registry
621 IOobject::NO_READ,
622 IOobject::AUTO_WRITE,
623 false
624 ),
625 s
626 );
627
628 Info<< "writing surfMesh:\n "
629 << smesh.searchableSurface::objectPath() << nl << endl;
630 smesh.searchableSurface::write();
631
632 Info<< "Written surface in = "
633 << timer.cpuTimeIncrement() << " s." << nl << endl;
634 }
635
636 mc.clean_all() ;
637
638 Info<< "End\n" << endl;
639
640 return 0;
641}
642
643
644// ************************************************************************* //
CGAL::Point_3< K > Point
CGAL::Triangle_3< K > Triangle
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
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
Random number generator.
Definition: Random.H:60
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
const labelList & surfaces() const
Return the surface indices.
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTimeCxx.H:54
A class for handling file names.
Definition: fileName.H:76
Identifies a surface patch/zone by name and index, with geometric type.
A triFace with additional (region) index.
Definition: labelledTri.H:60
Conversion functions between point (Foam::) and Point (CGAL::)
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
Implements a timeout mechanism via sigalarm.
Definition: timer.H:84
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
IOoject and searching on triSurface.
Triangulated surface description with patch information.
Definition: triSurface.H:79
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
rDeltaT ref()
const polyBoundaryMesh & patches
engineTime & runTime
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const dimensionedScalar c
Speed of light in a vacuum.
@ real
Definition: Roots.H:56
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
pointField vertices(const blockVertexList &bvl)
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
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
scalarField samples(nIntervals, Zero)
Random rndGen
Definition: createFields.H:23