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 |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2012-2016 OpenFOAM Foundation
9 Copyright (C) 2020-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
35\*---------------------------------------------------------------------------*/
36
37#include "argList.H"
38#include "Time.H"
39#include "searchableSurfaces.H"
41#include "triSurfaceMesh.H"
42#include "labelVector.H"
43
44#include "MarchingCubes.h"
45
46
47using namespace Foam;
48
49// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50
51// Main program:
52
53int main(int argc, char *argv[])
54{
55 argList::addNote
56 (
57 "Re-sample surfaces used in foamyHexMesh operation"
58 );
59 argList::addArgument("(nx ny nz)", "The resampling interval");
60 argList::addArgument("output", "The output triSurface/ file");
61
62 argList::noFunctionObjects(); // Never use function objects
63
64 #include "setRootCase.H"
65 #include "createTime.H"
66
67 const auto n = args.get<labelVector>(1);
68 const auto exportName = args.get<fileName>(2);
69
70 Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
71 << " writing re-sampled " << n << " to " << exportName
72 << nl << endl;
73
75
76 IOdictionary foamyHexMeshDict
77 (
79 (
80 "foamyHexMeshDict",
81 runTime.system(),
82 runTime,
83 IOobject::MUST_READ_IF_MODIFIED,
84 IOobject::NO_WRITE
85 )
86 );
87
88 // Define/load all geometry
89 searchableSurfaces allGeometry
90 (
92 (
93 "cvSearchableSurfaces",
94 runTime.constant(),
95 "triSurface",
96 runTime,
97 IOobject::MUST_READ,
98 IOobject::NO_WRITE
99 ),
100 foamyHexMeshDict.subDict("geometry"),
101 foamyHexMeshDict.getOrDefault("singleRegionName", true)
102 );
103
104 Info<< "Geometry read in = "
105 << timer.cpuTimeIncrement() << " s." << nl << endl;
106
107
108 Random rndGen(64293*Pstream::myProcNo());
109
110 conformationSurfaces geometryToConformTo
111 (
112 runTime,
113 rndGen,
114 allGeometry,
115 foamyHexMeshDict.subDict("surfaceConformation")
116 );
117
118 Info<< "Set up geometry in = "
119 << timer.cpuTimeIncrement() << " s." << nl << endl;
120
121
122
123 // Extend
124 treeBoundBox bb = geometryToConformTo.globalBounds();
125 {
126 const vector smallVec = 0.1*bb.span();
127 bb.min() -= smallVec;
128 bb.max() += smallVec;
129 }
130
131 Info<< "Meshing to bounding box " << bb << nl << endl;
132
133 const vector span(bb.span());
134 const vector d
135 (
136 span.x()/(n.x()-1),
137 span.y()/(n.y()-1),
138 span.z()/(n.z()-1)
139 );
140
141 MarchingCubes mc(span.x(), span.y(), span.z() ) ;
142 mc.set_resolution(n.x(), n.y(), n.z());
143 mc.init_all() ;
144
145
146 // Generate points
147 pointField points(mc.size_x()*mc.size_y()*mc.size_z());
148 label pointi = 0;
149
150 point pt;
151 for( int k = 0 ; k < mc.size_z() ; k++ )
152 {
153 pt.z() = bb.min().z() + k*d.z();
154 for( int j = 0 ; j < mc.size_y() ; j++ )
155 {
156 pt.y() = bb.min().y() + j*d.y();
157 for( int i = 0 ; i < mc.size_x() ; i++ )
158 {
159 pt.x() = bb.min().x() + i*d.x();
160 points[pointi++] = pt;
161 }
162 }
163 }
164
165 Info<< "Generated " << points.size() << " sampling points in = "
166 << timer.cpuTimeIncrement() << " s." << nl << endl;
167
168
169 // Compute data
170 const searchableSurfaces& geometry = geometryToConformTo.geometry();
171 const labelList& surfaces = geometryToConformTo.surfaces();
172
173 scalarField signedDist;
174 labelList nearestSurfaces;
175 searchableSurfacesQueries::signedDistance
176 (
177 geometry,
178 surfaces,
179 points,
180 scalarField(points.size(), sqr(GREAT)),
181 searchableSurface::OUTSIDE, // for non-closed surfaces treat as
182 // outside
183 nearestSurfaces,
184 signedDist
185 );
186
187 // Fill elements
188 pointi = 0;
189 for( int k = 0 ; k < mc.size_z() ; k++ )
190 {
191 for( int j = 0 ; j < mc.size_y() ; j++ )
192 {
193 for( int i = 0 ; i < mc.size_x() ; i++ )
194 {
195 mc.set_data(float(signedDist[pointi++]), i, j, k);
196 }
197 }
198 }
199
200 Info<< "Determined signed distance in = "
201 << timer.cpuTimeIncrement() << " s." << nl << endl;
202
203
204 mc.run() ;
205
206 Info<< "Constructed iso surface in = "
207 << timer.cpuTimeIncrement() << " s." << nl << endl;
208
209
210 mc.clean_temps() ;
211
212
213
214 // Write output file
215 if (mc.ntrigs() > 0)
216 {
217 Triangle* triangles = mc.triangles();
218 List<labelledTri> tris(mc.ntrigs());
219 forAll(tris, triI)
220 {
221 tris[triI] = labelledTri
222 (
223 triangles[triI].v1,
224 triangles[triI].v2,
225 triangles[triI].v3,
226 0 // region
227 );
228 }
229
230
231 Vertex* vertices = mc.vertices();
232 pointField points(mc.nverts());
233 forAll(points, pointi)
234 {
235 Vertex& v = vertices[pointi];
236 points[pointi] = point
237 (
238 bb.min().x() + v.x*span.x()/mc.size_x(),
239 bb.min().y() + v.y*span.y()/mc.size_y(),
240 bb.min().z() + v.z*span.z()/mc.size_z()
241 );
242 }
243
244
245 // Find correspondence to original surfaces
246 labelList regionOffsets(surfaces.size());
247 label nRegions = 0;
248 forAll(surfaces, i)
249 {
250 const wordList& regions = geometry[surfaces[i]].regions();
251 regionOffsets[i] = nRegions;
252 nRegions += regions.size();
253 }
254
255
257 nRegions = 0;
258 forAll(surfaces, i)
259 {
260 const wordList& regions = geometry[surfaces[i]].regions();
261
262 forAll(regions, regionI)
263 {
265 (
266 geometry[surfaces[i]].name() + "_" + regions[regionI],
267 nRegions,
268 "patch"
269 );
270 nRegions++;
271 }
272 }
273
274 triSurface s(tris, patches, points, true);
275
276 Info<< "Extracted triSurface in = "
277 << timer.cpuTimeIncrement() << " s." << nl << endl;
278
279
280 // Find out region on local surface of nearest point
281 {
282 List<pointIndexHit> hitInfo;
283 labelList hitSurfaces;
284 geometryToConformTo.findSurfaceNearest
285 (
286 s.faceCentres(),
287 scalarField(s.size(), sqr(GREAT)),
288 hitInfo,
289 hitSurfaces
290 );
291
292 // Get region
293 DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
294 DynamicList<label> surfIndices(hitSurfaces.size());
295
296 forAll(surfaces, surfI)
297 {
298 // Extract info on this surface
299 surfInfo.clear();
300 surfIndices.clear();
301 forAll(hitSurfaces, triI)
302 {
303 if (hitSurfaces[triI] == surfI)
304 {
305 surfInfo.append(hitInfo[triI]);
306 surfIndices.append(triI);
307 }
308 }
309
310 // Calculate sideness of these surface points
311 labelList region;
312 geometry[surfaces[surfI]].getRegion(surfInfo, region);
313
314 forAll(region, i)
315 {
316 label triI = surfIndices[i];
317 s[triI].region() = regionOffsets[surfI]+region[i];
318 }
319 }
320 }
321
322 Info<< "Re-patched surface in = "
323 << timer.cpuTimeIncrement() << " s." << nl << endl;
324
325 triSurfaceMesh smesh
326 (
328 (
329 exportName,
330 runTime.constant(), // instance
331 "triSurface",
332 runTime, // registry
333 IOobject::NO_READ,
334 IOobject::AUTO_WRITE,
335 false
336 ),
337 s
338 );
339
340 Info<< "writing surfMesh:\n "
341 << smesh.searchableSurface::objectPath() << nl << endl;
342 smesh.searchableSurface::write();
343
344 Info<< "Written surface in = "
345 << timer.cpuTimeIncrement() << " s." << nl << endl;
346 }
347
348 mc.clean_all() ;
349
350
351 Info<< "End\n" << endl;
352
353 return 0;
354}
355
356
357// ************************************************************************* //
CGAL::Triangle_3< K > Triangle
label k
graph_traits< Graph >::vertex_descriptor Vertex
Definition: SloanRenumber.C:74
label n
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
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
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
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
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
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...
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))
Namespace for OpenFOAM.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
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)
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
Random rndGen
Definition: createFields.H:23