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