PrimitivePatchProjectPoints.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) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2020-2021 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 Description
28  For every point on the patch find the closest face on the target side.
29  Return a target face label for each patch point
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "boolList.H"
34 #include "PointHit.H"
35 #include "objectHit.H"
36 #include "bandCompression.H"
37 
38 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
39 
40 template<class FaceList, class PointField>
41 template<class ToPatch>
44 (
45  const ToPatch& targetPatch,
46  const Field
47  <
49  >& projectionDirection,
50  const intersection::algorithm alg,
51  const intersection::direction dir
52 ) const
53 {
54  // The current patch is slave, i.e. it is being projected onto the target
55 
56  if (projectionDirection.size() != nPoints())
57  {
59  << "Projection direction field does not correspond to "
60  << "patch points." << endl
61  << "Size: " << projectionDirection.size()
62  << " Number of points: " << nPoints()
63  << abort(FatalError);
64  }
65 
66  const labelList& slavePointOrder = localPointOrder();
67 
68  const labelList& slaveMeshPoints = meshPoints();
69 
70  // Result
71  List<objectHit> result(nPoints());
72 
73  const labelListList& masterFaceFaces = targetPatch.faceFaces();
74 
75  const ToPatch& masterFaces = targetPatch;
76 
77  const Field<point_type>& masterPoints = targetPatch.points();
78 
79  // Estimate face centre of target side
80  Field<point_type> masterFaceCentres(targetPatch.size());
81 
82  forAll(masterFaceCentres, facei)
83  {
84  masterFaceCentres[facei] =
85  average(masterFaces[facei].points(masterPoints));
86  }
87 
88  // Algorithm:
89  // Loop through all points of the slave side. For every point find the
90  // radius for the current contact face. If the contact point falls inside
91  // the face and the radius is smaller than for all neighbouring faces,
92  // the contact is found. If not, visit the neighbour closest to the
93  // calculated contact point. If a single master face is visited more than
94  // twice, initiate n-squared search.
95 
96  label curFace = 0;
97  label nNSquaredSearches = 0;
98 
99  forAll(slavePointOrder, pointi)
100  {
101  // Pick up slave point and direction
102  const label curLocalPointLabel = slavePointOrder[pointi];
103 
104  const point_type& curPoint =
105  points_[slaveMeshPoints[curLocalPointLabel]];
106 
107  const point_type& curProjectionDir =
108  projectionDirection[curLocalPointLabel];
109 
110  bool closer;
111 
112  boolList visitedTargetFace(targetPatch.size(), false);
113  bool doNSquaredSearch = false;
114 
115  bool foundEligible = false;
116 
117  scalar sqrDistance = GREAT;
118 
119  // Force the full search for the first point to ensure good
120  // starting face
121  if (pointi == 0)
122  {
123  doNSquaredSearch = true;
124  }
125  else
126  {
127  do
128  {
129  closer = false;
130  doNSquaredSearch = false;
131 
132  // Calculate intersection with curFace
133  PointHit<point_type> curHit =
134  masterFaces[curFace].ray
135  (
136  curPoint,
137  curProjectionDir,
138  masterPoints,
139  alg,
140  dir
141  );
142 
143  visitedTargetFace[curFace] = true;
144 
145  if (curHit.hit())
146  {
147  result[curLocalPointLabel] = objectHit(true, curFace);
148 
149  break;
150  }
151  else
152  {
153  // If a new miss is eligible, it is closer than
154  // any previous eligible miss (due to surface walk)
155 
156  // Only grab the miss if it is eligible
157  if (curHit.eligibleMiss())
158  {
159  foundEligible = true;
160  result[curLocalPointLabel] = objectHit(false, curFace);
161  }
162 
163  // Find the next likely face for intersection
164 
165  // Calculate the miss point on the plane of the
166  // face. This is cooked (illogical!) for fastest
167  // surface walk.
168  //
169  point_type missPlanePoint =
170  curPoint + curProjectionDir*curHit.distance();
171 
172  const labelList& masterNbrs = masterFaceFaces[curFace];
173 
174  sqrDistance =
175  magSqr(missPlanePoint - masterFaceCentres[curFace]);
176 
177  forAll(masterNbrs, nbrI)
178  {
179  if
180  (
181  magSqr
182  (
183  missPlanePoint
184  - masterFaceCentres[masterNbrs[nbrI]]
185  )
186  <= sqrDistance
187  )
188  {
189  closer = true;
190  curFace = masterNbrs[nbrI];
191  }
192  }
193 
194  if (visitedTargetFace[curFace])
195  {
196  // This face has already been visited.
197  // Execute n-squared search
198  doNSquaredSearch = true;
199  break;
200  }
201  }
202 
203  DebugInfo << '.';
204  } while (closer);
205  }
206 
207  if
208  (
209  doNSquaredSearch || !foundEligible
210  )
211  {
212  nNSquaredSearches++;
213 
214  DebugInfo << "p " << curLocalPointLabel << ": ";
215 
216  result[curLocalPointLabel] = objectHit(false, -1);
217  scalar minDistance = GREAT;
218 
219  forAll(masterFaces, facei)
220  {
221  PointHit<point_type> curHit =
222  masterFaces[facei].ray
223  (
224  curPoint,
225  curProjectionDir,
226  masterPoints,
227  alg,
228  dir
229  );
230 
231  if (curHit.hit())
232  {
233  result[curLocalPointLabel] = objectHit(true, facei);
234  curFace = facei;
235 
236  break;
237  }
238  else if (curHit.eligibleMiss())
239  {
240  // Calculate min distance
241  scalar missDist =
242  Foam::mag(curHit.missPoint() - curPoint);
243 
244  if (missDist < minDistance)
245  {
246  minDistance = missDist;
247 
248  result[curLocalPointLabel] = objectHit(false, facei);
249  curFace = facei;
250  }
251  }
252  }
253 
254  DebugInfo << result[curLocalPointLabel] << nl;
255  }
256  else
257  {
258  DebugInfo << 'x';
259  }
260  }
261 
262  DebugInfo
263  << nl << "Executed " << nNSquaredSearches
264  << " n-squared searches out of total of "
265  << nPoints() << endl;
266 
267  return result;
268 }
269 
270 
271 template<class FaceList, class PointField>
272 template<class ToPatch>
275 (
276  const ToPatch& targetPatch,
277  const Field
278  <
280  >& projectionDirection,
281  const intersection::algorithm alg,
282  const intersection::direction dir
283 ) const
284 {
285  // The current patch is slave, i.e. it is being projected onto the target
286 
287  if (projectionDirection.size() != this->size())
288  {
290  << "Projection direction field does not correspond to patch faces."
291  << endl << "Size: " << projectionDirection.size()
292  << " Number of points: " << this->size()
293  << abort(FatalError);
294  }
295 
296  labelList slaveFaceOrder = bandCompression(faceFaces());
297 
298  // calculate master face centres
299  Field<point_type> masterFaceCentres(targetPatch.size());
300 
301  const labelListList& masterFaceFaces = targetPatch.faceFaces();
302 
303  const ToPatch& masterFaces = targetPatch;
304 
305  const typename ToPatch::PointFieldType& masterPoints = targetPatch.points();
306 
307  forAll(masterFaceCentres, facei)
308  {
309  masterFaceCentres[facei] =
310  masterFaces[facei].centre(masterPoints);
311  }
312 
313  // Result
314  List<objectHit> result(this->size());
315 
316  const PrimitivePatch<FaceList, PointField>& slaveFaces = *this;
317 
318  const PointField& slaveGlobalPoints = points();
319 
320  // Algorithm:
321  // Loop through all points of the slave side. For every point find the
322  // radius for the current contact face. If the contact point falls inside
323  // the face and the radius is smaller than for all neighbouring faces,
324  // the contact is found. If not, visit the neighbour closest to the
325  // calculated contact point. If a single master face is visited more than
326  // twice, initiate n-squared search.
327 
328  label curFace = 0;
329  label nNSquaredSearches = 0;
330 
331  forAll(slaveFaceOrder, facei)
332  {
333  // pick up slave point and direction
334  const label curLocalFaceLabel = slaveFaceOrder[facei];
335 
336  const point& curFaceCentre =
337  slaveFaces[curLocalFaceLabel].centre(slaveGlobalPoints);
338 
339  const vector& curProjectionDir =
340  projectionDirection[curLocalFaceLabel];
341 
342  bool closer;
343 
344  boolList visitedTargetFace(targetPatch.size(), false);
345  bool doNSquaredSearch = false;
346 
347  bool foundEligible = false;
348 
349  scalar sqrDistance = GREAT;
350 
351  // Force the full search for the first point to ensure good
352  // starting face
353  if (facei == 0)
354  {
355  doNSquaredSearch = true;
356  }
357  else
358  {
359  do
360  {
361  closer = false;
362  doNSquaredSearch = false;
363 
364  // Calculate intersection with curFace
365  PointHit<point_type> curHit =
366  masterFaces[curFace].ray
367  (
368  curFaceCentre,
369  curProjectionDir,
370  masterPoints,
371  alg,
372  dir
373  );
374 
375  visitedTargetFace[curFace] = true;
376 
377  if (curHit.hit())
378  {
379  result[curLocalFaceLabel] = objectHit(true, curFace);
380 
381  break;
382  }
383  else
384  {
385  // If a new miss is eligible, it is closer than
386  // any previous eligible miss (due to surface walk)
387 
388  // Only grab the miss if it is eligible
389  if (curHit.eligibleMiss())
390  {
391  foundEligible = true;
392  result[curLocalFaceLabel] = objectHit(false, curFace);
393  }
394 
395  // Find the next likely face for intersection
396 
397  // Calculate the miss point. This is
398  // cooked (illogical!) for fastest surface walk.
399  //
400  point_type missPlanePoint =
401  curFaceCentre + curProjectionDir*curHit.distance();
402 
403  sqrDistance =
404  magSqr(missPlanePoint - masterFaceCentres[curFace]);
405 
406  const labelList& masterNbrs = masterFaceFaces[curFace];
407 
408  forAll(masterNbrs, nbrI)
409  {
410  if
411  (
412  magSqr
413  (
414  missPlanePoint
415  - masterFaceCentres[masterNbrs[nbrI]]
416  )
417  <= sqrDistance
418  )
419  {
420  closer = true;
421  curFace = masterNbrs[nbrI];
422  }
423  }
424 
425  if (visitedTargetFace[curFace])
426  {
427  // This face has already been visited.
428  // Execute n-squared search
429  doNSquaredSearch = true;
430  break;
431  }
432  }
433 
434  DebugInfo << '.';
435  } while (closer);
436  }
437 
438  if (doNSquaredSearch || !foundEligible)
439  {
440  nNSquaredSearches++;
441 
442  DebugInfo << "p " << curLocalFaceLabel << ": ";
443 
444  result[curLocalFaceLabel] = objectHit(false, -1);
445  scalar minDistance = GREAT;
446 
447  forAll(masterFaces, facei)
448  {
449  PointHit<point_type> curHit =
450  masterFaces[facei].ray
451  (
452  curFaceCentre,
453  curProjectionDir,
454  masterPoints,
455  alg,
456  dir
457  );
458 
459  if (curHit.hit())
460  {
461  result[curLocalFaceLabel] = objectHit(true, facei);
462  curFace = facei;
463 
464  break;
465  }
466  else if (curHit.eligibleMiss())
467  {
468  // Calculate min distance
469  scalar missDist =
470  Foam::mag(curHit.missPoint() - curFaceCentre);
471 
472  if (missDist < minDistance)
473  {
474  minDistance = missDist;
475 
476  result[curLocalFaceLabel] = objectHit(false, facei);
477  curFace = facei;
478  }
479  }
480  }
481 
482  DebugInfo << result[curLocalFaceLabel] << nl;
483  }
484  else
485  {
486  DebugInfo << 'x';
487  }
488  }
489 
490  DebugInfo
491  << nl
492  << "Executed " << nNSquaredSearches
493  << " n-squared searches out of total of "
494  << this->size() << endl;
495 
496  return result;
497 }
498 
499 
500 // ************************************************************************* //
Foam::intersection::direction
direction
Definition: intersection.H:66
boolList.H
Foam::PointHit
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:53
Foam::bandCompression
labelList bandCompression(const labelListList &addressing)
Renumbers the addressing to reduce the band of the matrix.
Definition: bandCompression.C:44
Foam::PointHit::distance
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
objectHit.H
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::PointHit::eligibleMiss
bool eligibleMiss() const noexcept
Is this an eligible miss.
Definition: PointHit.H:127
Foam::fac::average
tmp< GeometricField< Type, faPatchField, areaMesh > > average(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Area-weighted average a edgeField creating a areaField.
Definition: facAverage.C:47
Foam::PointHit::missPoint
const point_type & missPoint() const
Return the miss point. Fatal if hit.
Definition: PointHit.H:158
Foam::PrimitivePatch::point_type
std::remove_reference< PointField >::type::value_type point_type
The point type.
Definition: PrimitivePatch.H:94
Foam::intersection::algorithm
algorithm
Definition: intersection.H:72
Foam::PrimitivePatch::projectFaceCentres
List< objectHit > projectFaceCentres(const ToPatch &targetPatch, const Field< point_type > &projectionDirection, const intersection::algorithm=intersection::FULL_RAY, const intersection::direction=intersection::VECTOR) const
Project vertices of patch onto another patch.
Foam::FatalError
error FatalError
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
bandCompression.H
The bandCompression function renumbers the addressing such that the band of the matrix is reduced....
Foam::PrimitivePatch::projectPoints
List< objectHit > projectPoints(const ToPatch &targetPatch, const Field< point_type > &projectionDirection, const intersection::algorithm=intersection::FULL_RAY, const intersection::direction=intersection::VECTOR) const
Project vertices of patch onto another patch.
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Vector< scalar >
Foam::List< Foam::objectHit >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
PointHit.H
Foam::objectHit
This class describes a combination of target object index and success flag. Behaves somewhat like std...
Definition: objectHit.H:50
Foam::PointHit::hit
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:79