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 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  const List<labelledTri>& tris = *this;
44  const pointField& points = this->points();
45 
46  Info<< nl << "# findLineAll did not hit its own face."
47  << nl << "# fi " << fi
48  << nl << "# start " << start
49  << nl << "# point " << p
50  << nl << "# end " << end
51  << nl << "# hitInfo " << hitInfo
52  << endl;
53 
54  meshTools::writeOBJ(Info, start);
57 
58  Info<< "l 1 2 3" << endl;
59 
60  meshTools::writeOBJ(Info, points[tris[fi][0]]);
61  meshTools::writeOBJ(Info, points[tris[fi][1]]);
62  meshTools::writeOBJ(Info, points[tris[fi][2]]);
63 
64  Info<< "f 4 5 6" << endl;
65 
66  forAll(hitInfo, hi)
67  {
68  label hfi = hitInfo[hi].index();
69 
70  meshTools::writeOBJ(Info, points[tris[hfi][0]]);
71  meshTools::writeOBJ(Info, points[tris[hfi][1]]);
72  meshTools::writeOBJ(Info, points[tris[hfi][2]]);
73 
74  Info<< "f "
75  << 3*hi + 7 << " "
76  << 3*hi + 8 << " "
77  << 3*hi + 9
78  << endl;
79  }
80 }
81 
82 
83 void Foam::triSurfaceMesh::processHit
84 (
85  scalar& internalCloseness,
86  scalar& externalCloseness,
87  const scalar internalToleranceCosAngle,
88  const scalar externalToleranceCosAngle,
89  const label fi,
90  const point& start,
91  const point& p,
92  const point& end,
93  const vector& normal,
94  const vectorField& normals,
95  const pointIndexHitList& hitInfo
96 ) const
97 {
98  if (hitInfo.size() < 1)
99  {
100  drawHitProblem(fi, start, p, end, hitInfo);
101  }
102  else if (hitInfo.size() == 1)
103  {
104  if (!hitInfo[0].hit())
105  {
106  }
107  else if (hitInfo[0].index() != fi)
108  {
109  drawHitProblem(fi, start, p, end, hitInfo);
110  }
111  }
112  else
113  {
114  label ownHiti = -1;
115 
116  forAll(hitInfo, hI)
117  {
118  // Find the hit on the triangle that launched the ray
119 
120  if (hitInfo[hI].index() == fi)
121  {
122  ownHiti = hI;
123  break;
124  }
125  }
126 
127  if (ownHiti < 0)
128  {
129  drawHitProblem(fi, start, p, end, hitInfo);
130  }
131  else if (ownHiti == 0)
132  {
133  // There are no internal hits, the first hit is the
134  // closest external hit
135 
136  if
137  (
138  (normal & normals[hitInfo[ownHiti + 1].index()])
139  < externalToleranceCosAngle
140  )
141  {
142  externalCloseness = min
143  (
144  externalCloseness,
145  mag(p - hitInfo[ownHiti + 1].hitPoint())
146  );
147  }
148  }
149  else if (ownHiti == hitInfo.size() - 1)
150  {
151  // There are no external hits, the last but one hit is
152  // the closest internal hit
153 
154  if
155  (
156  (normal & normals[hitInfo[ownHiti - 1].index()])
157  < internalToleranceCosAngle
158  )
159  {
160  internalCloseness = min
161  (
162  internalCloseness,
163  mag(p - hitInfo[ownHiti - 1].hitPoint())
164  );
165  }
166  }
167  else
168  {
169  if
170  (
171  (normal & normals[hitInfo[ownHiti + 1].index()])
172  < externalToleranceCosAngle
173  )
174  {
175  externalCloseness = min
176  (
177  externalCloseness,
178  mag(p - hitInfo[ownHiti + 1].hitPoint())
179  );
180  }
181 
182  if
183  (
184  (normal & normals[hitInfo[ownHiti - 1].index()])
185  < internalToleranceCosAngle
186  )
187  {
188  internalCloseness = min
189  (
190  internalCloseness,
191  mag(p - hitInfo[ownHiti - 1].hitPoint())
192  );
193  }
194  }
195  }
196 }
197 
198 
201 (
202  const scalar internalAngleTolerance,
203  const scalar externalAngleTolerance
204 ) const
205 {
206  const scalar internalToleranceCosAngle
207  (
208  cos(degToRad(180 - internalAngleTolerance))
209  );
210 
211  const scalar externalToleranceCosAngle
212  (
213  cos(degToRad(180 - externalAngleTolerance))
214  );
215 
216  const Time& runTime = objectRegistry::time();
217 
218  // Prepare start and end points for intersection tests
219 
220  const vectorField& normals = faceNormals();
221 
222  const scalar span = bounds().mag();
223 
224  const pointField start(faceCentres() - span*normals);
225  const pointField end(faceCentres() + span*normals);
226  const pointField& faceCentres = this->faceCentres();
227 
228  List<List<pointIndexHit>> allHitinfo;
229 
230  // Find all intersections (in order)
231  findLineAll(start, end, allHitinfo);
232 
233  scalarField internalCloseness(start.size(), great);
234  scalarField externalCloseness(start.size(), great);
235 
236  forAll(allHitinfo, fi)
237  {
238  const List<pointIndexHit>& hitInfo = allHitinfo[fi];
239 
240  processHit
241  (
242  internalCloseness[fi],
243  externalCloseness[fi],
244  internalToleranceCosAngle,
245  externalToleranceCosAngle,
246  fi,
247  start[fi],
248  faceCentres[fi],
249  end[fi],
250  normals[fi],
251  normals,
252  hitInfo
253  );
254  }
255 
257  (
259  (
261  (
262  IOobject
263  (
264  objectRegistry::name() + ".internalCloseness",
265  runTime.constant(),
266  "triSurface",
267  runTime
268  ),
269  *this,
270  dimLength,
271  internalCloseness
272  )
273  ),
274 
276  (
278  (
279  IOobject
280  (
281  objectRegistry::name() + ".externalCloseness",
282  runTime.constant(),
283  "triSurface",
284  runTime
285  ),
286  *this,
287  dimLength,
288  externalCloseness
289  )
290  )
291  );
292 }
293 
294 
297 (
298  const scalar internalAngleTolerance,
299  const scalar externalAngleTolerance
300 ) const
301 {
302  const scalar internalToleranceCosAngle
303  (
304  cos(degToRad(180 - internalAngleTolerance))
305  );
306 
307  const scalar externalToleranceCosAngle
308  (
309  cos(degToRad(180 - externalAngleTolerance))
310  );
311 
312  const Time& runTime = objectRegistry::time();
313 
314  // Prepare start and end points for intersection tests
315 
316  const pointField& points = this->points();
317  const labelList& meshPoints = this->meshPoints();
318  const pointField& faceCentres = this->faceCentres();
319  const vectorField& normals = this->faceNormals();
320  const labelListList& pointFaces = this->pointFaces();
321 
322  const scalar span = bounds().mag();
323 
324  label nPointFaces = 0;
325  forAll(pointFaces, pfi)
326  {
327  nPointFaces += pointFaces[pfi].size();
328  }
329 
330  pointField facePoints(nPointFaces);
331  pointField start(nPointFaces);
332  pointField end(nPointFaces);
333 
334  label i = 0;
335  forAll(points, pi)
336  {
337  forAll(pointFaces[pi], pfi)
338  {
339  const label fi = pointFaces[pi][pfi];
340 
341  facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]);
342  const vector& n = normals[fi];
343 
344  start[i] = facePoints[i] - span*n;
345  end[i] = facePoints[i] + span*n;
346 
347  i++;
348  }
349  }
350 
351  List<pointIndexHitList> allHitinfo;
352 
353  // Find all intersections (in order)
354  findLineAll(start, end, allHitinfo);
355 
356  scalarField internalCloseness(points.size(), great);
357  scalarField externalCloseness(points.size(), great);
358 
359  i = 0;
360  forAll(points, pi)
361  {
362  forAll(pointFaces[pi], pfi)
363  {
364  const label fi = pointFaces[pi][pfi];
365  const pointIndexHitList& hitInfo = allHitinfo[i];
366 
367  processHit
368  (
369  internalCloseness[pi],
370  externalCloseness[pi],
371  internalToleranceCosAngle,
372  externalToleranceCosAngle,
373  fi,
374  start[i],
375  facePoints[i],
376  end[i],
377  normals[fi],
378  normals,
379  hitInfo
380  );
381 
382  i++;
383  }
384  }
385 
387  (
389  (
391  (
392  IOobject
393  (
394  objectRegistry::name() + ".internalPointCloseness",
395  runTime.constant(),
396  "triSurface",
397  runTime
398  ),
399  *this,
400  dimLength,
401  internalCloseness
402  )
403  ),
404 
406  (
408  (
409  IOobject
410  (
411  objectRegistry::name() + ".externalPointCloseness",
412  runTime.constant(),
413  "triSurface",
414  runTime
415  ),
416  *this,
417  dimLength,
418  externalCloseness
419  )
420  )
421  );
422 }
423 
424 
425 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
const word & name() const
Return name.
Definition: IOobject.H:297
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:110
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:60
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:163
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const Field< point > & faceCentres() const
Return face centres for patch.
const Field< point > & faceNormals() const
Return face normals for patch.
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:68
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.
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)
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:90
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit >> &) const
Get all intersections in order from start to end.
const labelListList & pointFaces() const
Return point-face addressing.
Pair< tmp< triSurfacePointScalarField > > extractPointCloseness(const scalar internalAngleTolerance=80, const scalar externalAngleTolerance=10) const
Return a pair of triSurfaceScalarPointFields representing the.
static const char nl
Definition: Ostream.H:265
const Time & time() const
Return time.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
vector point
Point is a vector.
Definition: point.H:41
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Pair< tmp< triSurfaceScalarField > > extractCloseness(const scalar internalAngleTolerance=80, const scalar externalAngleTolerance=10) const
Return a pair of triSurfaceScalarFields representing the.
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:92
Foam::DimensionedField< scalar, triSurfaceGeoMesh > triSurfaceScalarField