FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/boolList.H>
31 #include <OpenFOAM/PointHit_.H>
32 #include <OpenFOAM/objectHit.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  const PointField& slaveGlobalPoints = points();
339 
340  // Algorithm:
341  // Loop through all points of the slave side. For every point find the
342  // radius for the current contact face. If the contact point falls inside
343  // the face and the radius is smaller than for all neighbouring faces,
344  // the contact is found. If not, visit the neighbour closest to the
345  // calculated contact point. If a single master face is visited more than
346  // twice, initiate n-squared search.
347 
348  label curFace = 0;
349  label nNSquaredSearches = 0;
350 
351  forAll (slaveFaceOrder, faceI)
352  {
353  // pick up slave point and direction
354  const label curLocalFaceLabel = slaveFaceOrder[faceI];
355 
356  const point& curFaceCentre =
357  slaveFaces[curLocalFaceLabel].centre(slaveGlobalPoints);
358 
359  const vector& curProjectionDir =
360  projectionDirection[curLocalFaceLabel];
361 
362  bool closer;
363 
364  boolList visitedTargetFace(targetPatch.size(), false);
365  bool doNSquaredSearch = false;
366 
367  bool foundEligible = false;
368 
369  scalar sqrDistance = GREAT;
370 
371  // Force the full search for the first point to ensure good
372  // starting face
373  if (faceI == 0)
374  {
375  doNSquaredSearch = true;
376  }
377  else
378  {
379  do
380  {
381  closer = false;
382  doNSquaredSearch = false;
383 
384  // Calculate intersection with curFace
385  PointHit<PointType> curHit =
386  masterFaces[curFace].ray
387  (
388  curFaceCentre,
389  curProjectionDir,
390  masterPoints,
391  alg,
392  dir
393  );
394 
395  visitedTargetFace[curFace] = true;
396 
397  if (curHit.hit())
398  {
399  result[curLocalFaceLabel] = objectHit(true, curFace);
400 
401  break;
402  }
403  else
404  {
405  // If a new miss is eligible, it is closer than
406  // any previous eligible miss (due to surface walk)
407 
408  // Only grab the miss if it is eligible
409  if (curHit.eligibleMiss())
410  {
411  foundEligible = true;
412  result[curLocalFaceLabel] = objectHit(false, curFace);
413  }
414 
415  // Find the next likely face for intersection
416 
417  // Calculate the miss point. This is
418  // cooked (illogical!) for fastest surface walk.
419  //
420  PointType missPlanePoint =
421  curFaceCentre + curProjectionDir*curHit.distance();
422 
423  sqrDistance =
424  magSqr(missPlanePoint - masterFaceCentres[curFace]);
425 
426  const labelList& masterNbrs = masterFaceFaces[curFace];
427 
428  forAll (masterNbrs, nbrI)
429  {
430  if
431  (
432  magSqr
433  (
434  missPlanePoint
435  - masterFaceCentres[masterNbrs[nbrI]]
436  )
437  <= sqrDistance
438  )
439  {
440  closer = true;
441  curFace = masterNbrs[nbrI];
442  }
443  }
444 
445  if (visitedTargetFace[curFace])
446  {
447  // This face has already been visited.
448  // Execute n-squared search
449  doNSquaredSearch = true;
450  break;
451  }
452  }
453 
454  if (debug) Info<< ".";
455  } while (closer);
456  }
457 
458  if (doNSquaredSearch || !foundEligible)
459  {
460  nNSquaredSearches++;
461 
462  if (debug)
463  {
464  Info<< "p " << curLocalFaceLabel << ": ";
465  }
466 
467  result[curLocalFaceLabel] = objectHit(false, -1);
468  scalar minDistance = GREAT;
469 
470  forAll (masterFaces, faceI)
471  {
472  PointHit<PointType> curHit =
473  masterFaces[faceI].ray
474  (
475  curFaceCentre,
476  curProjectionDir,
477  masterPoints,
478  alg,
479  dir
480  );
481 
482  if (curHit.hit())
483  {
484  result[curLocalFaceLabel] = objectHit(true, faceI);
485  curFace = faceI;
486 
487  break;
488  }
489  else if (curHit.eligibleMiss())
490  {
491  // Calculate min distance
492  scalar missDist =
493  Foam::mag(curHit.missPoint() - curFaceCentre);
494 
495  if (missDist < minDistance)
496  {
497  minDistance = missDist;
498 
499  result[curLocalFaceLabel] = objectHit(false, faceI);
500  curFace = faceI;
501  }
502  }
503  }
504 
505  if (debug)
506  {
507  Info << result[curLocalFaceLabel] << nl;
508  }
509  }
510  else
511  {
512  if (debug) Info<< "x";
513  }
514  }
515 
516  if (debug)
517  {
518  Info<< nl << "Executed " << nNSquaredSearches
519  << " n-squared searches out of total of "
520  << this->size() << endl;
521  }
522 
523  return result;
524 }
525 
526 
527 // ************************ vim: set sw=4 sts=4 et: ************************ //