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 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 
35 \*---------------------------------------------------------------------------*/
36 
37 #include "argList.H"
38 #include "Time.H"
39 #include "searchableSurfaces.H"
40 #include "conformationSurfaces.H"
41 #include "triSurfaceMesh.H"
42 #include "labelVector.H"
43 
44 #include "MarchingCubes.h"
45 
46 
47 using namespace Foam;
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 // Main program:
52 
53 int main(int argc, char *argv[])
54 {
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 labelVector n(args.get<labelVector>(1));
68  const fileName exportName = args[2];
69 
70  Info<< "Reading surfaces as specified in the foamyHexMeshDict and"
71  << " writing re-sampled " << n << " to " << exportName
72  << nl << endl;
73 
74  cpuTime timer;
75 
76  IOdictionary foamyHexMeshDict
77  (
78  IOobject
79  (
80  "foamyHexMeshDict",
81  runTime.system(),
82  runTime,
85  )
86  );
87 
88  // Define/load all geometry
89  searchableSurfaces allGeometry
90  (
91  IOobject
92  (
93  "cvSearchableSurfaces",
94  runTime.constant(),
95  "triSurface",
96  runTime,
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 
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;
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  {
264  patches[nRegions] = geometricSurfacePatch
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  (
327  IOobject
328  (
329  exportName,
330  runTime.constant(), // instance
331  "triSurface",
332  runTime, // registry
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 // ************************************************************************* //
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
Foam::Vector::x
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
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::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
labelVector.H
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Vertex
graph_traits< Graph >::vertex_descriptor Vertex
Definition: SloanRenumber.C:75
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
Foam::boundBox::max
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
Foam::argList::get
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:251
Foam::conformationSurfaces
Definition: conformationSurfaces.H:55
Foam::Vector::z
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
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
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::argList::noFunctionObjects
static void noFunctionObjects(bool addWithOption=false)
Remove '-noFunctionObjects' option and ignore any occurrences.
Definition: argList.C:467
Foam::Field< vector >
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
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::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
Time.H
setRootCase.H
Foam::Vector::y
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
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< label >
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
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
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
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::IOobject::MUST_READ
Definition: IOobject.H:120