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