extractCloseness.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) 2018-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 \*---------------------------------------------------------------------------*/
25 
26 #include "triSurfaceMesh.H"
27 #include "triSurfaceFields.H"
28 #include "meshTools.H"
29 #include "Time.H"
30 #include "unitConversion.H"
31 
32 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
33 
34 void Foam::triSurfaceMesh::drawHitProblem
35 (
36  const label fi,
37  const point& start,
38  const point& p,
39  const point& end,
40  const pointIndexHitList& hitInfo
41 ) const
42 {
43  if (debug)
44  {
45  const List<labelledTri>& tris = *this;
46  const pointField& points = this->points();
47 
48  Info<< nl << "# findLineAll did not hit its own face."
49  << nl << "# fi " << fi
50  << nl << "# start " << start
51  << nl << "# point " << p
52  << nl << "# end " << end
53  << nl << "# hitInfo " << hitInfo
54  << endl;
55 
56  meshTools::writeOBJ(Info, start);
59 
60  Info<< "l 1 2 3" << endl;
61 
62  meshTools::writeOBJ(Info, points[tris[fi][0]]);
63  meshTools::writeOBJ(Info, points[tris[fi][1]]);
64  meshTools::writeOBJ(Info, points[tris[fi][2]]);
65 
66  Info<< "f 4 5 6" << endl;
67 
68  forAll(hitInfo, hi)
69  {
70  label hfi = hitInfo[hi].index();
71 
72  meshTools::writeOBJ(Info, points[tris[hfi][0]]);
73  meshTools::writeOBJ(Info, points[tris[hfi][1]]);
74  meshTools::writeOBJ(Info, points[tris[hfi][2]]);
75 
76  Info<< "f "
77  << 3*hi + 7 << " "
78  << 3*hi + 8 << " "
79  << 3*hi + 9
80  << endl;
81  }
82  }
83 }
84 
85 
86 void Foam::triSurfaceMesh::processHit
87 (
88  scalar& internalCloseness,
89  scalar& externalCloseness,
90  const scalar internalToleranceCosAngle,
91  const scalar externalToleranceCosAngle,
92  const label fi,
93  const point& start,
94  const point& p,
95  const point& end,
96  const vector& normal,
97  const vectorField& normals,
98  const pointIndexHitList& hitInfo
99 ) const
100 {
101  if (hitInfo.size() < 1)
102  {
103  drawHitProblem(fi, start, p, end, hitInfo);
104  }
105  else if (hitInfo.size() == 1)
106  {
107  if (!hitInfo[0].hit())
108  {
109  }
110  else if (hitInfo[0].index() != fi)
111  {
112  drawHitProblem(fi, start, p, end, hitInfo);
113  }
114  }
115  else
116  {
117  label ownHiti = -1;
118 
119  forAll(hitInfo, hI)
120  {
121  // Find the hit on the triangle that launched the ray
122 
123  if (hitInfo[hI].index() == fi)
124  {
125  ownHiti = hI;
126  break;
127  }
128  }
129 
130  if (ownHiti < 0)
131  {
132  drawHitProblem(fi, start, p, end, hitInfo);
133  }
134  else if (ownHiti == 0)
135  {
136  // There are no internal hits, the first hit is the
137  // closest external hit
138 
139  if
140  (
141  (normal & normals[hitInfo[ownHiti + 1].index()])
142  < externalToleranceCosAngle
143  )
144  {
145  externalCloseness = min
146  (
147  externalCloseness,
148  mag(p - hitInfo[ownHiti + 1].hitPoint())
149  );
150  }
151  }
152  else if (ownHiti == hitInfo.size() - 1)
153  {
154  // There are no external hits, the last but one hit is
155  // the closest internal hit
156 
157  if
158  (
159  (normal & normals[hitInfo[ownHiti - 1].index()])
160  < internalToleranceCosAngle
161  )
162  {
163  internalCloseness = min
164  (
165  internalCloseness,
166  mag(p - hitInfo[ownHiti - 1].hitPoint())
167  );
168  }
169  }
170  else
171  {
172  if
173  (
174  (normal & normals[hitInfo[ownHiti + 1].index()])
175  < externalToleranceCosAngle
176  )
177  {
178  externalCloseness = min
179  (
180  externalCloseness,
181  mag(p - hitInfo[ownHiti + 1].hitPoint())
182  );
183  }
184 
185  if
186  (
187  (normal & normals[hitInfo[ownHiti - 1].index()])
188  < internalToleranceCosAngle
189  )
190  {
191  internalCloseness = min
192  (
193  internalCloseness,
194  mag(p - hitInfo[ownHiti - 1].hitPoint())
195  );
196  }
197  }
198  }
199 }
200 
201 
204 (
205  const scalar internalAngleTolerance,
206  const scalar externalAngleTolerance
207 ) const
208 {
209  const scalar internalToleranceCosAngle
210  (
211  cos(degToRad(180 - internalAngleTolerance))
212  );
213 
214  const scalar externalToleranceCosAngle
215  (
216  cos(degToRad(180 - externalAngleTolerance))
217  );
218 
219  const Time& runTime = objectRegistry::time();
220 
221  // Prepare start and end points for intersection tests
222 
223  const vectorField& normals = faceNormals();
224 
225  const scalar span = bounds().mag();
226 
227  const pointField start(faceCentres() - span*normals);
228  const pointField end(faceCentres() + span*normals);
229  const pointField& faceCentres = this->faceCentres();
230 
231  List<List<pointIndexHit>> allHitinfo;
232 
233  // Find all intersections (in order)
234  findLineAll(start, end, allHitinfo);
235 
236  scalarField internalCloseness(start.size(), great);
237  scalarField externalCloseness(start.size(), great);
238 
239  forAll(allHitinfo, fi)
240  {
241  const List<pointIndexHit>& hitInfo = allHitinfo[fi];
242 
243  processHit
244  (
245  internalCloseness[fi],
246  externalCloseness[fi],
247  internalToleranceCosAngle,
248  externalToleranceCosAngle,
249  fi,
250  start[fi],
251  faceCentres[fi],
252  end[fi],
253  normals[fi],
254  normals,
255  hitInfo
256  );
257  }
258 
260  (
262  (
264  (
265  IOobject
266  (
267  objectRegistry::name() + ".internalCloseness",
268  runTime.constant(),
270  runTime
271  ),
272  *this,
273  dimLength,
274  internalCloseness
275  )
276  ),
277 
279  (
281  (
282  IOobject
283  (
284  objectRegistry::name() + ".externalCloseness",
285  runTime.constant(),
287  runTime
288  ),
289  *this,
290  dimLength,
291  externalCloseness
292  )
293  )
294  );
295 }
296 
297 
300 (
301  const scalar internalAngleTolerance,
302  const scalar externalAngleTolerance
303 ) const
304 {
305  const scalar internalToleranceCosAngle
306  (
307  cos(degToRad(180 - internalAngleTolerance))
308  );
309 
310  const scalar externalToleranceCosAngle
311  (
312  cos(degToRad(180 - externalAngleTolerance))
313  );
314 
315  const Time& runTime = objectRegistry::time();
316 
317  // Prepare start and end points for intersection tests
318 
319  const pointField& points = this->points();
320  const labelList& meshPoints = this->meshPoints();
321  const pointField& faceCentres = this->faceCentres();
322  const vectorField& normals = this->faceNormals();
323  const labelListList& pointFaces = this->pointFaces();
324 
325  const scalar span = bounds().mag();
326 
327  label nPointFaces = 0;
328  forAll(pointFaces, pfi)
329  {
330  nPointFaces += pointFaces[pfi].size();
331  }
332 
333  pointField facePoints(nPointFaces);
334  pointField start(nPointFaces);
335  pointField end(nPointFaces);
336 
337  label i = 0;
338  forAll(points, pi)
339  {
340  forAll(pointFaces[pi], pfi)
341  {
342  const label fi = pointFaces[pi][pfi];
343 
344  facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]);
345  const vector& n = normals[fi];
346 
347  start[i] = facePoints[i] - span*n;
348  end[i] = facePoints[i] + span*n;
349 
350  i++;
351  }
352  }
353 
354  List<pointIndexHitList> allHitinfo;
355 
356  // Find all intersections (in order)
357  findLineAll(start, end, allHitinfo);
358 
359  scalarField internalCloseness(points.size(), great);
360  scalarField externalCloseness(points.size(), great);
361 
362  i = 0;
363  forAll(points, pi)
364  {
365  forAll(pointFaces[pi], pfi)
366  {
367  const label fi = pointFaces[pi][pfi];
368  const pointIndexHitList& hitInfo = allHitinfo[i];
369 
370  processHit
371  (
372  internalCloseness[pi],
373  externalCloseness[pi],
374  internalToleranceCosAngle,
375  externalToleranceCosAngle,
376  fi,
377  start[i],
378  facePoints[i],
379  end[i],
380  normals[fi],
381  normals,
382  hitInfo
383  );
384 
385  i++;
386  }
387  }
388 
390  (
392  (
394  (
395  IOobject
396  (
397  objectRegistry::name() + ".internalPointCloseness",
398  runTime.constant(),
400  runTime
401  ),
402  *this,
403  dimLength,
404  internalCloseness
405  )
406  ),
407 
409  (
411  (
412  IOobject
413  (
414  objectRegistry::name() + ".externalPointCloseness",
415  runTime.constant(),
417  runTime
418  ),
419  *this,
420  dimLength,
421  externalCloseness
422  )
423  )
424  );
425 }
426 
427 
428 // ************************************************************************* //
Pair< tmp< triSurfacePointScalarField > > extractPointCloseness(const scalar internalAngleTolerance=80, const scalar externalAngleTolerance=80) const
Return a pair of triSurfaceScalarPointFields representing the.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
Definition: IOobject.H:315
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
const Field< PointType > & faceCentres() const
Return face centres for patch.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Fields for triSurface.
virtual tmp< pointField > points() const
Get the points that define the surface.
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:69
List< pointIndexHit > pointIndexHitList
List of pointIndexHits and associated functions.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
const dimensionSet dimLength
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dimensionedScalar cos(const dimensionedScalar &ds)
static const word & geometryDir()
Return the geometry directory name.
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:90
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit >> &) const
Get all intersections in order from start to end.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
const Field< PointType > & faceNormals() const
Return face normals for patch.
Pair< tmp< triSurfaceScalarField > > extractCloseness(const scalar internalAngleTolerance=80, const scalar externalAngleTolerance=80) const
Return a pair of triSurfaceScalarFields representing the.
const labelListList & pointFaces() const
Return point-face addressing.
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
vector point
Point is a vector.
Definition: point.H:41
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
const boundBox & bounds() const
Return const reference to boundBox.
Foam::DimensionedField< scalar, triSurfacePointGeoMesh > triSurfacePointScalarField
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
Foam::DimensionedField< scalar, triSurfaceGeoMesh > triSurfaceScalarField