All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
foamyHexMeshSurfaceSimplify_non_octree.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-2021 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 
32 \*---------------------------------------------------------------------------*/
33 
34 #include "argList.H"
35 #include "Time.H"
36 #include "searchableSurfaces.H"
37 #include "conformationSurfaces.H"
38 #include "triSurfaceMesh.H"
39 
40 #include "MarchingCubes.h"
41 
42 
43 using namespace Foam;
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
47 // Main program:
48 
49 int main(int argc, char *argv[])
50 {
52  (
53  "Re-sample surfaces used in foamyHexMesh operation"
54  );
55  // argList::validArgs.append("inputFile");
56  argList::validArgs.append("(nx ny nz)");
57  argList::validArgs.append("outputName");
58 
59  #include "setRootCase.H"
60  #include "createTime.H"
62 
63  const Vector<label> n(IStringStream(args.args()[1])());
64  const fileName exportName = args.args()[2];
65 
66  Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
67  << " writing re-sampled " << n << " to " << exportName
68  << nl << endl;
69 
70  cpuTime timer;
71 
72  IOdictionary foamyHexMeshDict
73  (
74  IOobject
75  (
76  "foamyHexMeshDict",
77  runTime.system(),
78  runTime,
81  )
82  );
83 
84  // Define/load all geometry
85  searchableSurfaces allGeometry
86  (
87  IOobject
88  (
89  "cvSearchableSurfaces",
90  runTime.constant(),
92  runTime,
95  ),
96  foamyHexMeshDict.subDict("geometry"),
97  foamyHexMeshDict.lookupOrDefault("singleRegionName", true)
98  );
99 
100  Info<< "Geometry read in = "
101  << timer.cpuTimeIncrement() << " s." << nl << endl;
102 
103 
105 
106  conformationSurfaces geometryToConformTo
107  (
108  runTime,
109  rndGen,
110  allGeometry,
111  foamyHexMeshDict.subDict("surfaceConformation")
112  );
113 
114  Info<< "Set up geometry in = "
115  << timer.cpuTimeIncrement() << " s." << nl << endl;
116 
117 
118 
119  // Extend
120  treeBoundBox bb = geometryToConformTo.globalBounds();
121  {
122  const vector smallVec = 0.1*bb.span();
123  bb.min() -= smallVec;
124  bb.max() += smallVec;
125  }
126 
127  Info<< "Meshing to bounding box " << bb << nl << endl;
128 
129  const vector span(bb.span());
130  const vector d
131  (
132  span.x()/(n.x()-1),
133  span.y()/(n.y()-1),
134  span.z()/(n.z()-1)
135  );
136 
137  MarchingCubes mc(span.x(), span.y(), span.z() ) ;
138  mc.set_resolution(n.x(), n.y(), n.z());
139  mc.init_all() ;
140 
141 
142  // Generate points
143  pointField points(mc.size_x()*mc.size_y()*mc.size_z());
144  label pointi = 0;
145 
146  point pt;
147  for( int k = 0 ; k < mc.size_z() ; k++ )
148  {
149  pt.z() = bb.min().z() + k*d.z();
150  for( int j = 0 ; j < mc.size_y() ; j++ )
151  {
152  pt.y() = bb.min().y() + j*d.y();
153  for( int i = 0 ; i < mc.size_x() ; i++ )
154  {
155  pt.x() = bb.min().x() + i*d.x();
156  points[pointi++] = pt;
157  }
158  }
159  }
160 
161  Info<< "Generated " << points.size() << " sampling points in = "
162  << timer.cpuTimeIncrement() << " s." << nl << endl;
163 
164 
165  // Compute data
166  const searchableSurfaces& geometry = geometryToConformTo.geometry();
167  const labelList& surfaces = geometryToConformTo.surfaces();
168 
169  scalarField signedDist;
170  labelList nearestSurfaces;
172  (
173  geometry,
174  surfaces,
175  points,
176  scalarField(points.size(), sqr(great)),
177  searchableSurface::OUTSIDE, // for non-closed surfaces treat as
178  // outside
179  nearestSurfaces,
180  signedDist
181  );
182 
183  // Fill elements
184  pointi = 0;
185  for( int k = 0 ; k < mc.size_z() ; k++ )
186  {
187  for( int j = 0 ; j < mc.size_y() ; j++ )
188  {
189  for( int i = 0 ; i < mc.size_x() ; i++ )
190  {
191  mc.set_data(float(signedDist[pointi++]), i, j, k);
192  }
193  }
194  }
195 
196  Info<< "Determined signed distance in = "
197  << timer.cpuTimeIncrement() << " s." << nl << endl;
198 
199 
200  mc.run() ;
201 
202  Info<< "Constructed iso surface in = "
203  << timer.cpuTimeIncrement() << " s." << nl << endl;
204 
205 
206  mc.clean_temps() ;
207 
208 
209 
210  // Write output file
211  if (mc.ntrigs() > 0)
212  {
213  Triangle* triangles = mc.triangles();
214  List<labelledTri> tris(mc.ntrigs());
215  forAll(tris, triI)
216  {
217  tris[triI] = labelledTri
218  (
219  triangles[triI].v1,
220  triangles[triI].v2,
221  triangles[triI].v3,
222  0 // region
223  );
224  }
225 
226 
227  Vertex* vertices = mc.vertices();
228  pointField points(mc.nverts());
229  forAll(points, pointi)
230  {
231  Vertex& v = vertices[pointi];
232  points[pointi] = point
233  (
234  bb.min().x() + v.x*span.x()/mc.size_x(),
235  bb.min().y() + v.y*span.y()/mc.size_y(),
236  bb.min().z() + v.z*span.z()/mc.size_z()
237  );
238  }
239 
240 
241  // Find correspondence to original surfaces
242  labelList regionOffsets(surfaces.size());
243  label nRegions = 0;
244  forAll(surfaces, i)
245  {
246  const wordList& regions = geometry[surfaces[i]].regions();
247  regionOffsets[i] = nRegions;
248  nRegions += regions.size();
249  }
250 
251 
253  nRegions = 0;
254  forAll(surfaces, i)
255  {
256  const wordList& regions = geometry[surfaces[i]].regions();
257 
258  forAll(regions, regionI)
259  {
260  patches[nRegions] = geometricSurfacePatch
261  (
262  "patch",
263  geometry[surfaces[i]].name() + "_" + regions[regionI],
264  nRegions
265  );
266  nRegions++;
267  }
268  }
269 
270  triSurface s(tris, patches, points, true);
271 
272  Info<< "Extracted triSurface in = "
273  << timer.cpuTimeIncrement() << " s." << nl << endl;
274 
275 
276  // Find out region on local surface of nearest point
277  {
278  List<pointIndexHit> hitInfo;
279  labelList hitSurfaces;
280  geometryToConformTo.findSurfaceNearest
281  (
282  s.faceCentres(),
283  scalarField(s.size(), sqr(great)),
284  hitInfo,
285  hitSurfaces
286  );
287 
288  // Get region
289  DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
290  DynamicList<label> surfIndices(hitSurfaces.size());
291 
292  forAll(surfaces, surfI)
293  {
294  // Extract info on this surface
295  surfInfo.clear();
296  surfIndices.clear();
297  forAll(hitSurfaces, triI)
298  {
299  if (hitSurfaces[triI] == surfI)
300  {
301  surfInfo.append(hitInfo[triI]);
302  surfIndices.append(triI);
303  }
304  }
305 
306  // Calculate sideness of these surface points
307  labelList region;
308  geometry[surfaces[surfI]].getRegion(surfInfo, region);
309 
310  forAll(region, i)
311  {
312  label triI = surfIndices[i];
313  s[triI].region() = regionOffsets[surfI]+region[i];
314  }
315  }
316  }
317 
318  Info<< "Re-patched surface in = "
319  << timer.cpuTimeIncrement() << " s." << nl << endl;
320 
321  triSurfaceMesh smesh
322  (
323  IOobject
324  (
325  exportName,
326  runTime.constant(),
328  runTime,
331  false
332  ),
333  s
334  );
335 
336  Info<< "writing surfMesh:\n "
337  << smesh.searchableSurface::objectPath() << nl << endl;
338  smesh.searchableSurface::write();
339 
340  Info<< "Written surface in = "
341  << timer.cpuTimeIncrement() << " s." << nl << endl;
342  }
343 
344  mc.clean_all() ;
345 
346 
347  Info<< "End\n" << endl;
348 
349  return 0;
350 }
351 
352 
353 // ************************************************************************* //
graph_traits< Graph >::vertex_descriptor Vertex
Definition: SloanRenumber.C:72
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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:79
void off()
Switch the function objects off.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
const Cmpt & z() const
Definition: VectorI.H:87
label k
Boltzmann constant.
patches[0]
const Cmpt & y() const
Definition: VectorI.H:81
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:53
IOoject and searching on triSurface.
Random rndGen(label(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
static const word & geometryDir()
Return the geometry directory name.
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:123
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:113
static const char nl
Definition: Ostream.H:260
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
vector point
Point is a vector.
Definition: point.H:41
const functionObjectList & functionObjects() const
Return the list of function objects.
Definition: Time.H:410
Input from memory buffer stream.
Definition: IStringStream.H:49
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
label n
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:159
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:72