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-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
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 Application
28  foamyHexMeshSurfaceSimplify
29 
30 Description
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"
41 #include "conformationSurfaces.H"
42 #include "triSurfaceMesh.H"
43 
44 #include "opt_octree.h"
45 #include "cube.h"
46 
47 using namespace Foam;
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 class pointConversion
52 {
53  const vector scale_;
54 
55  const vector offset_;
56 
57 public:
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
87 class 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 
113  (
114  geometry,
115  surfaces,
116  samples,
117  scalarField(1, GREAT),
119  nearestSurfaces,
120  distance
121  );
122 
123  return distance[0];
124  }
125 
126 
127 public:
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 
367 int main(int argc, char *argv[])
368 {
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 fileName exportName = args[1];
382 
383  Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
384  << " writing a re-sampled surface to " << exportName
385  << nl << endl;
386 
387  cpuTime timer;
388 
389  IOdictionary foamyHexMeshDict
390  (
391  IOobject
392  (
393  "foamyHexMeshDict",
394  runTime.system(),
395  runTime,
398  )
399  );
400 
401  // Define/load all geometry
402  searchableSurfaces allGeometry
403  (
404  IOobject
405  (
406  "cvSearchableSurfaces",
407  runTime.constant(),
408  "triSurface",
409  runTime,
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 
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 =
460  cmptDivide
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  {
551  patches[nRegions] = geometricSurfacePatch
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  (
615  IOobject
616  (
617  exportName,
618  runTime.constant(), // instance
619  "triSurface",
620  runTime, // registry
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 // ************************************************************************* //
Foam::IOobject::NO_WRITE
Definition: IOobject.H:130
Foam::IOdictionary
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:54
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::Random
Random number generator.
Definition: Random.H:59
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:104
pointConversion
Conversion functions between point (Foam::) and Point (CGAL::)
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:129
Triangle
CGAL::Triangle_3< K > Triangle
Definition: CGALIndexedPolyhedron.H:57
Foam::fileName
A class for handling file names.
Definition: fileName.H:69
s
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))
Definition: gmvOutputSpray.H:25
Foam::cmptMultiply
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:55
Foam::geometricSurfacePatch
The geometricSurfacePatch is like patchIdentifier but for surfaces. Holds type, name and index.
Definition: geometricSurfacePatch.H:52
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Foam::argList::addNote
static void addNote(const string &note)
Add extra notes for the usage information.
Definition: argList.C:426
Foam::triSurfaceMesh
IOoject and searching on triSurface.
Definition: triSurfaceMesh.H:106
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
ref
rDeltaT ref()
Foam::conformationSurfaces
Definition: conformationSurfaces.H:55
Foam::boundBox::min
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::argList::addArgument
static void addArgument(const string &argName, const string &usage="")
Append a (mandatory) argument to validArgs.
Definition: argList.C:315
Foam::boundBox::span
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
searchableSurfaces.H
Foam::argList::noFunctionObjects
static void noFunctionObjects(bool addWithOption=false)
Remove '-noFunctionObjects' option and ignore any occurrences.
Definition: argList.C:467
Foam::Field< scalar >
Foam::triSurface
Triangulated surface description with patch information.
Definition: triSurface.H:76
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:474
argList.H
Foam::searchableSurfacesQueries::signedDistance
static void signedDistance(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, const volumeType illegalHandling, labelList &nearestSurfaces, scalarField &distance)
Find signed distance to nearest surface. Outside is positive.
Definition: searchableSurfacesQueries.C:599
conformationSurfaces.H
Foam::conformationSurfaces::geometry
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
Definition: conformationSurfacesI.H:30
samples
scalarField samples(nIntervals, Zero)
Foam::roots::real
Definition: Roots.H:56
Foam::cpuTimeCxx
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTimeCxx.H:53
Foam::vertices
pointField vertices(const blockVertexList &bvl)
Definition: blockVertexList.H:49
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::conformationSurfaces::surfaces
const labelList & surfaces() const
Return the surface indices.
Definition: conformationSurfacesI.H:49
Foam::distance
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Foam::cmptDivide
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
Time.H
setRootCase.H
Foam::TimePaths::system
const word & system() const
Return system name.
Definition: TimePathsI.H:94
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:464
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::Vector< scalar >
Foam::List< label >
Foam::searchableSurfaces
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
Definition: searchableSurfaces.H:92
Foam::timer
Implements a timeout mechanism via sigalarm.
Definition: timer.H:83
points
const pointField & points
Definition: gmvOutputHeader.H:1
createTime.H
Foam::IOobject::MUST_READ_IF_MODIFIED
Definition: IOobject.H:121
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
Foam::labelledTri
Triangle with additional region number.
Definition: labelledTri.H:57
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Point
CGAL::Point_3< K > Point
Definition: CGALIndexedPolyhedron.H:53
rndGen
Random rndGen
Definition: createFields.H:23
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::TimePaths::constant
const word & constant() const
Return constant name.
Definition: TimePathsI.H:88
Foam::IOobject::NO_READ
Definition: IOobject.H:123
args
Foam::argList args(argc, argv)
triSurfaceMesh.H
Foam::volumeType::OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
Foam::IOobject::MUST_READ
Definition: IOobject.H:120