foamyHexMeshSurfaceSimplify.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  foamyHexMeshSurfaceSimplify
26 
27 Description
28  Simplifies surfaces by resampling.
29 
30  Uses Thomas Lewiner's topology preserving MarchingCubes.
31  (http://zeus.mat.puc-rio.br/tomlew/tomlew_uk.php)
32 
33 \*---------------------------------------------------------------------------*/
34 
35 #include "argList.H"
36 #include "Time.H"
37 #include "searchableSurfaces.H"
38 #include "conformationSurfaces.H"
39 #include "triSurfaceMesh.H"
40 
41 #include "opt_octree.h"
42 #include "cube.h"
43 
44 using namespace Foam;
45 
46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 
48 class pointConversion
49 {
50  const vector scale_;
51 
52  const vector offset_;
53 
54 public:
55 
56  //- Construct from components
58  (
59  const vector scale,
60  const vector offset
61  )
62  :
63  scale_(scale),
64  offset_(offset)
65  {}
66 
67  inline Point toLocal(const Foam::point& pt) const
68  {
69  Foam::point p = cmptMultiply(scale_, (pt + offset_));
70  return Point(p.x(), p.y(), p.z());
71  }
72 
73  inline Foam::point toGlobal(const Point& pt) const
74  {
75  point p(pt.x(), pt.y(), pt.z());
76  return cmptDivide(p, scale_) - offset_;
77  }
78 };
79 
80 
81 
82 
83 // For use in Fast-Dual Octree from Thomas Lewiner
84 class distanceCalc
85 :
86  public ::data_access
87 {
88 
89  const Level min_level_;
90 
91  const conformationSurfaces& geometryToConformTo_;
92 
93  const pointConversion& converter_;
94 
95 
96  // Private Member Functions
97 
98  scalar signedDistance(const Foam::point& pt) const
99  {
100  const searchableSurfaces& geometry = geometryToConformTo_.geometry();
101  const labelList& surfaces = geometryToConformTo_.surfaces();
102 
103  static labelList nearestSurfaces;
104  static scalarField distance;
105 
106  static pointField samples(1);
107  samples[0] = pt;
108 
110  (
111  geometry,
112  surfaces,
113  samples,
114  scalarField(1, great),
116  nearestSurfaces,
117  distance
118  );
119 
120  return distance[0];
121  }
122 
123 
124 public:
125 
126  // Constructors
127 
128  //- Construct from components
129  distanceCalc
130  (
131  Level max_level_,
132  real iso_val_,
133  Level min_level,
134  const conformationSurfaces& geometryToConformTo,
135  const pointConversion& converter
136  )
137  :
138  data_access(max_level_,iso_val_),
139  min_level_(min_level),
140  geometryToConformTo_(geometryToConformTo),
141  converter_(converter)
142  {}
143 
144 
145  //- Destructor
146  virtual ~distanceCalc()
147  {}
148 
149 
150  // Member Functions
151 
152  //- Test function
153  virtual bool need_refine( const Cube &c )
154  {
155  int l = c.lv() ;
156 
157  if ( l >= _max_level ) return false;
158  if ( l < min_level_ ) return true;
159 
160  treeBoundBox bb
161  (
162  converter_.toGlobal
163  (
164  Point
165  (
166  c.xmin(),
167  c.ymin(),
168  c.zmin()
169  )
170  ),
171  converter_.toGlobal
172  (
173  Point
174  (
175  c.xmax(),
176  c.ymax(),
177  c.zmax()
178  )
179  )
180  );
181 
182  const searchableSurfaces& geometry =
183  geometryToConformTo_.geometry();
184  const labelList& surfaces =
185  geometryToConformTo_.surfaces();
186 
187 
188  //- Uniform refinement around surface
189  {
190  forAll(surfaces, i)
191  {
192  if (geometry[surfaces[i]].overlaps(bb))
193  {
194  return true;
195  }
196  }
197  return false;
198  }
199 
200 
202  // scalar ccDist = signedDistance(bb.midpoint());
203  // scalar ccVal = ccDist - _iso_val;
204  // if (mag(ccVal) < small)
205  //{
206  // return true;
207  //}
208  // const pointField points(bb.points());
209  // forAll(points, pointi)
210  //{
211  // scalar pointVal = signedDistance(points[pointi]) - _iso_val;
212  // if (ccVal*pointVal < 0)
213  // {
214  // return true;
215  // }
216  //}
217  // return false;
218 
219 
223  // const pointField points(bb.points());
224  // const edgeList& edges = treeBoundBox::edges;
225  // pointField start(edges.size());
226  // pointField end(edges.size());
227  // forAll(edges, i)
228  //{
229  // start[i] = points[edges[i][0]];
230  // end[i] = points[edges[i][1]];
231  //}
232  // Foam::List<Foam::List<pointIndexHit>> hitInfo;
233  // labelListList hitSurfaces;
234  // searchableSurfacesQueries::findAllIntersections
235  //(
236  // geometry,
237  // surfaces,
238  // start,
239  // end,
240  // hitSurfaces,
241  // hitInfo
242  //);
243  //
245  // label nInt = 0;
246  // forAll(hitSurfaces, edgeI)
247  //{
248  // nInt += hitSurfaces[edgeI].size();
249  //}
250  //
251  // if (nInt == 0)
252  //{
253  // // No surface intersected. See if there is one inside
254  // forAll(surfaces, i)
255  // {
256  // if (geometry[surfaces[i]].overlaps(bb))
257  // {
258  // return true;
259  // }
260  // }
261  // return false;
262  //}
263  //
265  // label baseSurfI = -1;
266  // forAll(hitSurfaces, edgeI)
267  //{
268  // const labelList& hSurfs = hitSurfaces[edgeI];
269  // forAll(hSurfs, i)
270  // {
271  // if (baseSurfI == -1)
272  // {
273  // baseSurfI = hSurfs[i];
274  // }
275  // else if (baseSurfI != hSurfs[i])
276  // {
277  // // Multiple surfaces
278  // return true;
279  // }
280  // }
281  //}
282  //
284  // DynamicList<pointIndexHit> baseInfo(nInt);
285  // forAll(hitInfo, edgeI)
286  //{
287  // const Foam::List<pointIndexHit>& hits = hitInfo[edgeI];
288  // forAll(hits, i)
289  // {
290  // (void)hits[i].hitPoint();
291  // baseInfo.append(hits[i]);
292  // }
293  //}
294  // vectorField normals;
295  // geometry[surfaces[baseSurfI]].getNormal(baseInfo, normals);
296  // for (label i = 1; i < normals.size(); ++i)
297  //{
298  // if ((normals[0] & normals[i]) < 0.9)
299  // {
300  // return true;
301  // }
302  //}
303  // labelList regions;
304  // geometry[surfaces[baseSurfI]].getRegion(baseInfo, regions);
305  // for (label i = 1; i < regions.size(); ++i)
306  //{
307  // if (regions[0] != regions[i])
308  // {
309  // return true;
310  // }
311  //}
312  // return false;
313 
314 
315 
316  // samples[0] = point(c.xmin(), c.ymin(), c.zmin());
317  // samples[1] = point(c.xmax(), c.ymin(), c.zmin());
318  // samples[2] = point(c.xmax(), c.ymax(), c.zmin());
319  // samples[3] = point(c.xmin(), c.ymax(), c.zmin());
320  //
321  // samples[4] = point(c.xmin(), c.ymin(), c.zmax());
322  // samples[5] = point(c.xmax(), c.ymin(), c.zmax());
323  // samples[6] = point(c.xmax(), c.ymax(), c.zmax());
324  // samples[7] = point(c.xmin(), c.ymax(), c.zmax());
325 
326  // scalarField nearestDistSqr(8, great);
327  //
328  // Foam::List<pointIndexHit> nearestInfo;
329  // surf_.findNearest(samples, nearestDistSqr, nearestInfo);
330  // vectorField normals;
331  // surf_.getNormal(nearestInfo, normals);
332  //
333  // for (label i = 1; i < normals.size(); ++i)
334  //{
335  // if ((normals[0] & normals[i]) < 0.5)
336  // {
337  // return true;
338  // }
339  //}
340  // return false;
341 
343  // const labelList elems(surf_.tree().findBox(bb));
344  //
345  // if (elems.size() > 1)
346  //{
347  // return true;
348  //}
349  // else
350  //{
351  // return false;
352  //}
353  }
354 
355  //- Data function
356  virtual real value_at( const Cube &c )
357  {
358  return signedDistance(converter_.toGlobal(c)) - _iso_val;
359  }
360 };
361 
362 
363 // Main program:
364 
365 int main(int argc, char *argv[])
366 {
368  (
369  "Re-sample surfaces used in foamyHexMesh operation"
370  );
371  argList::validArgs.append("outputName");
372 
373  #include "setRootCase.H"
374  #include "createTime.H"
376 
377  const fileName exportName = args.args()[1];
378 
379  Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
380  << " writing a re-sampled surface to " << exportName
381  << nl << endl;
382 
383  cpuTime timer;
384 
385  IOdictionary foamyHexMeshDict
386  (
387  IOobject
388  (
389  "foamyHexMeshDict",
390  runTime.system(),
391  runTime,
394  )
395  );
396 
397  // Define/load all geometry
398  searchableSurfaces allGeometry
399  (
400  IOobject
401  (
402  "cvSearchableSurfaces",
403  runTime.constant(),
404  "triSurface",
405  runTime,
408  ),
409  foamyHexMeshDict.subDict("geometry"),
410  foamyHexMeshDict.lookupOrDefault("singleRegionName", true)
411  );
412 
413  Info<< "Geometry read in = "
414  << timer.cpuTimeIncrement() << " s." << nl << endl;
415 
416 
418 
419  conformationSurfaces geometryToConformTo
420  (
421  runTime,
422  rndGen,
423  allGeometry,
424  foamyHexMeshDict.subDict("surfaceConformation")
425  );
426 
427  Info<< "Set up geometry in = "
428  << timer.cpuTimeIncrement() << " s." << nl << endl;
429 
430 
431  const searchableSurfaces& geometry = geometryToConformTo.geometry();
432  const labelList& surfaces = geometryToConformTo.surfaces();
433 
434 
435  const label minLevel = 2;
436 
437  // The max cube size follows from the minLevel and the default cube size
438  // (1)
439  const scalar maxSize = 1.0 / (1 << minLevel);
440  const scalar halfMaxSize = 0.5*maxSize;
441 
442 
443  // Scale the geometry to fit within
444  // halfMaxSize .. 1-halfMaxSize
445 
446  scalar wantedRange = 1.0-maxSize;
447 
448  const treeBoundBox bb = geometryToConformTo.globalBounds();
449 
450  const vector scale = cmptDivide
451  (
452  point(wantedRange, wantedRange, wantedRange),
453  bb.span()
454  );
455  const vector offset =
456  cmptDivide
457  (
458  point(halfMaxSize, halfMaxSize, halfMaxSize),
459  scale
460  )
461  -bb.min();
462 
463 
464  const pointConversion converter(scale, offset);
465 
466 
467  // Marching cubes
468 
469  OptOctree octree;
470 
471  distanceCalc ref
472  (
473  8, // maxLevel
474  0.0, // distance
475  minLevel, // minLevel
476  geometryToConformTo,
477  converter
478  );
479 
480  octree.refine(&ref);
481  octree.set_impl(&ref);
482 
483  Info<< "Calculated octree in = "
484  << timer.cpuTimeIncrement() << " s." << nl << endl;
485 
486  MarchingCubes& mc = octree.mc();
487 
488  mc.clean_all() ;
489  octree.build_isosurface(&ref) ;
490 
491  Info<< "Constructed iso surface of distance in = "
492  << timer.cpuTimeIncrement() << " s." << nl << endl;
493 
494  // Write output file
495  if (mc.ntrigs() > 0)
496  {
497  Triangle* triangles = mc.triangles();
498  label nTris = mc.ntrigs();
499  Foam::DynamicList<labelledTri> tris(mc.ntrigs());
500  for (label triI = 0; triI < nTris; ++triI)
501  {
502  const Triangle& t = triangles[triI];
503  if (t.v1 != t.v2 && t.v1 != t.v3 && t.v2 != t.v3)
504  {
505  tris.append
506  (
508  (
509  triangles[triI].v1,
510  triangles[triI].v2,
511  triangles[triI].v3,
512  0 // region
513  )
514  );
515  }
516  }
517 
518 
519  Point* vertices = mc.vertices();
520  pointField points(mc.nverts());
521  forAll(points, pointi)
522  {
523  const Point& v = vertices[pointi];
524  points[pointi] = converter.toGlobal(v);
525  }
526 
527 
528  // Find correspondence to original surfaces
529  labelList regionOffsets(surfaces.size());
530  label nRegions = 0;
531  forAll(surfaces, i)
532  {
533  const wordList& regions = geometry[surfaces[i]].regions();
534  regionOffsets[i] = nRegions;
535  nRegions += regions.size();
536  }
537 
538 
540  nRegions = 0;
541  forAll(surfaces, i)
542  {
543  const wordList& regions = geometry[surfaces[i]].regions();
544 
545  forAll(regions, regionI)
546  {
547  patches[nRegions] = geometricSurfacePatch
548  (
549  "patch",
550  geometry[surfaces[i]].name() + "_" + regions[regionI],
551  nRegions
552  );
553  nRegions++;
554  }
555  }
556 
557  triSurface s(tris, patches, points, true);
558  tris.clearStorage();
559 
560  Info<< "Extracted triSurface in = "
561  << timer.cpuTimeIncrement() << " s." << nl << endl;
562 
563 
564  // Find out region on local surface of nearest point
565  {
567  labelList hitSurfaces;
568  geometryToConformTo.findSurfaceNearest
569  (
570  s.faceCentres(),
571  scalarField(s.size(), sqr(great)),
572  hitInfo,
573  hitSurfaces
574  );
575 
576  // Get region
577  DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
578  DynamicList<label> surfIndices(hitSurfaces.size());
579 
580  forAll(surfaces, surfI)
581  {
582  // Extract info on this surface
583  surfInfo.clear();
584  surfIndices.clear();
585  forAll(hitSurfaces, triI)
586  {
587  if (hitSurfaces[triI] == surfI)
588  {
589  surfInfo.append(hitInfo[triI]);
590  surfIndices.append(triI);
591  }
592  }
593 
594  // Calculate sideness of these surface points
595  labelList region;
596  geometry[surfaces[surfI]].getRegion(surfInfo, region);
597 
598  forAll(region, i)
599  {
600  label triI = surfIndices[i];
601  s[triI].region() = regionOffsets[surfI]+region[i];
602  }
603  }
604  }
605 
606  Info<< "Re-patched surface in = "
607  << timer.cpuTimeIncrement() << " s." << nl << endl;
608 
609  triSurfaceMesh smesh
610  (
611  IOobject
612  (
613  exportName,
614  runTime.constant(), // instance
615  "triSurface",
616  runTime, // registry
619  false
620  ),
621  s
622  );
623 
624  Info<< "writing surfMesh:\n "
625  << smesh.searchableSurface::objectPath() << nl << endl;
626  smesh.searchableSurface::write();
627 
628  Info<< "Written surface in = "
629  << timer.cpuTimeIncrement() << " s." << nl << endl;
630  }
631 
632  mc.clean_all() ;
633 
634  Info<< "End\n" << endl;
635 
636  return 0;
637 }
638 
639 
640 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A class for handling file names.
Definition: fileName.H:69
void off()
Switch the function objects off.
const labelList & surfaces() const
Return the surface indices.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:427
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
const Cmpt & z() const
Definition: VectorI.H:87
patches[0]
const Cmpt & y() const
Definition: VectorI.H:81
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
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.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
IOoject and searching on triSurface.
Random rndGen(label(0))
Conversion functions between point (Foam::) and Point (CGAL::)
scalarField samples(nIntervals, 0)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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 pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:74
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
Container for searchableSurfaces.
Triangle with additional region number.
Definition: labelledTri.H:57
const Cmpt & x() const
Definition: VectorI.H:75
Random number generator.
Definition: Random.H:57
const word & system() const
Return system name.
Definition: TimePaths.H:114
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:265
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
vector point
Point is a vector.
Definition: point.H:41
const functionObjectList & functionObjects() const
Return the list of function objects.
Definition: Time.H:410
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:84
The geometricSurfacePatch is like patchIdentifier but for surfaces. Holds type, name and index...
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
messageStream Info
rDeltaT ref()
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTime.H:52
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
static void addNote(const string &)
Add extra notes for the usage information.
Definition: argList.C:151
volScalarField & p
Triangulated surface description with patch information.
Definition: triSurface.H:66
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
const stringList & args() const
Return arguments.
Definition: argListI.H:66