FreeFOAM The Cross-Platform CFD Toolkit
searchableSurfacesQueries.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 \*----------------------------------------------------------------------------*/
25 
27 #include <OpenFOAM/ListOps.H>
28 #include <OpenFOAM/OFstream.H>
29 #include <meshTools/meshTools.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 Foam::pointIndexHit Foam::searchableSurfacesQueries::tempFindNearest
38 (
39  const searchableSurface& surf,
40  const point& pt,
41  const scalar initDistSqr
42 )
43 {
44  pointField onePoint(1, pt);
45  scalarField oneDist(1, initDistSqr);
46  List<pointIndexHit> oneHit(1);
47  surf.findNearest(onePoint, oneDist, oneHit);
48  return oneHit[0];
49 }
50 
51 
52 // Calculate sum of distance to surfaces.
53 Foam::scalar Foam::searchableSurfacesQueries::sumDistSqr
54 (
55  const PtrList<searchableSurface>& allSurfaces,
56  const labelList& surfacesToTest,
57  const scalar initDistSqr,
58  const point& pt
59 )
60 {
61  scalar sum = 0;
62 
63  forAll(surfacesToTest, testI)
64  {
65  label surfI = surfacesToTest[testI];
66 
67  pointIndexHit hit
68  (
69  tempFindNearest(allSurfaces[surfI], pt, initDistSqr)
70  );
71 
72  // Note: make it fall over if not hit.
73  sum += magSqr(hit.hitPoint()-pt);
74  }
75  return sum;
76 }
77 
78 
79 // Reflects the point furthest away around the triangle centre by a factor fac.
80 // (triangle centre is the average of all points but the ihi. pSum is running
81 // sum of all points)
82 Foam::scalar Foam::searchableSurfacesQueries::tryMorphTet
83 (
84  const PtrList<searchableSurface>& allSurfaces,
85  const labelList& surfacesToTest,
86  const scalar initDistSqr,
87  List<vector>& p,
88  List<scalar>& y,
89  vector& pSum,
90  const label ihi,
91  const scalar fac
92 )
93 {
94  scalar fac1 = (1.0-fac)/vector::nComponents;
95  scalar fac2 = fac1-fac;
96 
97  vector ptry = pSum*fac1-p[ihi]*fac2;
98 
99  scalar ytry = sumDistSqr(allSurfaces, surfacesToTest, initDistSqr, ptry);
100 
101  if (ytry < y[ihi])
102  {
103  y[ihi] = ytry;
104  pSum += ptry - p[ihi];
105  p[ihi] = ptry;
106  }
107  return ytry;
108 }
109 
110 
111 bool Foam::searchableSurfacesQueries::morphTet
112 (
113  const PtrList<searchableSurface>& allSurfaces,
114  const labelList& surfacesToTest,
115  const scalar initDistSqr,
116  const scalar convergenceDistSqr,
117  const label maxIter,
118  List<vector>& p,
119  List<scalar>& y
120 )
121 {
122  vector pSum = sum(p);
123 
124  autoPtr<OFstream> str;
125  label vertI = 0;
126  if (debug)
127  {
128  wordList names(surfacesToTest.size());
129  forAll(surfacesToTest, i)
130  {
131  names[i] = allSurfaces[surfacesToTest[i]].name();
132  }
133  Pout<< "searchableSurfacesQueries::morphTet : intersection of "
134  << names << " starting from points:" << p << endl;
135  str.reset(new OFstream("track.obj"));
136  meshTools::writeOBJ(str(), p[0]);
137  vertI++;
138  }
139 
140  for (label iter = 0; iter < maxIter; iter++)
141  {
142  // Get the indices of lowest, highest and second-highest values.
143  label ilo, ihi, inhi;
144  {
145  labelList sortedIndices;
146  sortedOrder(y, sortedIndices);
147  ilo = sortedIndices[0];
148  ihi = sortedIndices[sortedIndices.size()-1];
149  inhi = sortedIndices[sortedIndices.size()-2];
150  }
151 
152  if (debug)
153  {
154  Pout<< "Iteration:" << iter
155  << " lowest:" << y[ilo] << " highest:" << y[ihi]
156  << " points:" << p << endl;
157 
158  meshTools::writeOBJ(str(), p[ilo]);
159  vertI++;
160  str()<< "l " << vertI-1 << ' ' << vertI << nl;
161  }
162 
163  if (y[ihi] < convergenceDistSqr)
164  {
165  // Get point on 0th surface.
166  Swap(p[0], p[ilo]);
167  Swap(y[0], y[ilo]);
168  return true;
169  }
170 
171  // Reflection: point furthest away gets reflected.
172  scalar ytry = tryMorphTet
173  (
174  allSurfaces,
175  surfacesToTest,
176  10*y[ihi], // search box.
177  p,
178  y,
179  pSum,
180  ihi,
181  -1.0
182  );
183 
184  if (ytry <= y[ilo])
185  {
186  // If in right direction (y lower) expand by two.
187  ytry = tryMorphTet
188  (
189  allSurfaces,
190  surfacesToTest,
191  10*y[ihi],
192  p,
193  y,
194  pSum,
195  ihi,
196  2.0
197  );
198  }
199  else if (ytry >= y[inhi])
200  {
201  // If inside tet try contraction.
202 
203  scalar ysave = y[ihi];
204 
205  ytry = tryMorphTet
206  (
207  allSurfaces,
208  surfacesToTest,
209  10*y[ihi],
210  p,
211  y,
212  pSum,
213  ihi,
214  0.5
215  );
216 
217  if (ytry >= ysave)
218  {
219  // Contract around lowest point.
220  forAll(p, i)
221  {
222  if (i != ilo)
223  {
224  p[i] = 0.5*(p[i] + p[ilo]);
225  y[i] = sumDistSqr
226  (
227  allSurfaces,
228  surfacesToTest,
229  y[ihi],
230  p[i]
231  );
232  }
233  }
234  pSum = sum(p);
235  }
236  }
237  }
238 
239  if (debug)
240  {
241  meshTools::writeOBJ(str(), p[0]);
242  vertI++;
243  str()<< "l " << vertI-1 << ' ' << vertI << nl;
244  }
245 
246  // Failure to converge. Return best guess so far.
247  label ilo = findMin(y);
248  Swap(p[0], p[ilo]);
249  Swap(y[0], y[ilo]);
250  return false;
251 }
252 
253 
255 //void Foam::searchableSurfacesQueries::findAllIntersections
256 //(
257 // const searchableSurface& s,
258 // const pointField& start,
259 // const pointField& end,
260 // const vectorField& smallVec,
261 // List<List<pointIndexHit> >& surfaceHitInfo
262 //)
263 //{
264 // surfaceHitInfo.setSize(start.size());
265 //
266 // // Current start point of vector
267 // pointField p0(start);
268 //
269 // List<pointIndexHit> intersectInfo(start.size());
270 //
271 // // For test whether finished doing vector.
272 // const vectorField dirVec(end-start);
273 // const scalarField magSqrDirVec(magSqr(dirVec));
274 //
275 // while (true)
276 // {
277 // // Find first intersection. Synced.
278 // s.findLine(p0, end, intersectInfo);
279 //
280 // label nHits = 0;
281 //
282 // forAll(intersectInfo, i)
283 // {
284 // if (intersectInfo[i].hit())
285 // {
286 // nHits++;
287 //
288 // label sz = surfaceHitInfo[i].size();
289 // surfaceHitInfo[i].setSize(sz+1);
290 // surfaceHitInfo[i][sz] = intersectInfo[i];
291 //
292 // p0[i] = intersectInfo[i].hitPoint() + smallVec[i];
293 //
294 // // If beyond endpoint set to endpoint so as not to pick up
295 // // any intersections. Could instead just filter out hits.
296 // if (((p0[i]-start[i])&dirVec[i]) > magSqrDirVec[i])
297 // {
298 // p0[i] = end[i];
299 // }
300 // }
301 // else
302 // {
303 // // Set to endpoint to stop intersection test. See above.
304 // p0[i] = end[i];
305 // }
306 // }
307 //
308 // // returnReduce(nHits) ?
309 // if (nHits == 0)
310 // {
311 // break;
312 // }
313 // }
314 //}
315 
316 
317 // Given current set of hits (allSurfaces, allInfo) merge in those coming from
318 // surface surfI.
319 void Foam::searchableSurfacesQueries::mergeHits
320 (
321  const point& start,
322  const scalar mergeDist,
323 
324  const label testI, // index of surface
325  const List<pointIndexHit>& surfHits, // hits on surface
326 
327  labelList& allSurfaces,
328  List<pointIndexHit>& allInfo,
329  scalarList& allDistSqr
330 )
331 {
332  // Precalculate distances
333  scalarList surfDistSqr(surfHits.size());
334  forAll(surfHits, i)
335  {
336  surfDistSqr[i] = magSqr(surfHits[i].hitPoint()-start);
337  }
338 
339  forAll(surfDistSqr, i)
340  {
341  label index = findLower(allDistSqr, surfDistSqr[i]);
342 
343  // Check if equal to lower.
344  if
345  (
346  index >= 0
347  && (mag(allDistSqr[index]-surfDistSqr[i]) < mergeDist)
348  )
349  {
350  // Same. Do not count.
351  //Pout<< "point:" << surfHits[i].hitPoint()
352  // << " considered same as:" << allInfo[index].hitPoint()
353  // << " within tol:" << mergeDist
354  // << endl;
355  }
356  else
357  {
358  // Check if equal to higher
359  label next = index+1;
360  if
361  (
362  next < allDistSqr.size()
363  && (mag(allDistSqr[next]-surfDistSqr[i]) < mergeDist)
364  )
365  {
366  //Pout<< "point:" << surfHits[i].hitPoint()
367  // << " considered same as:" << allInfo[next].hitPoint()
368  // << " within tol:" << mergeDist
369  // << endl;
370  }
371  else
372  {
373  // Insert after index
374  label sz = allSurfaces.size();
375  allSurfaces.setSize(sz+1);
376  allInfo.setSize(allSurfaces.size());
377  allDistSqr.setSize(allSurfaces.size());
378  // Make space.
379  for (label j = sz-1; j > index; --j)
380  {
381  allSurfaces[j+1] = allSurfaces[j];
382  allInfo[j+1] = allInfo[j];
383  allDistSqr[j+1] = allDistSqr[j];
384  }
385  // Insert new value
386  allSurfaces[index+1] = testI;
387  allInfo[index+1] = surfHits[i];
388  allDistSqr[index+1] = surfDistSqr[i];
389  }
390  }
391  }
392 }
393 
394 
395 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
396 
397 // Find any intersection
399 (
400  const PtrList<searchableSurface>& allSurfaces,
401  const labelList& surfacesToTest,
402  const pointField& start,
403  const pointField& end,
404  labelList& hitSurfaces,
405  List<pointIndexHit>& hitInfo
406 )
407 {
408  hitSurfaces.setSize(start.size());
409  hitSurfaces = -1;
410  hitInfo.setSize(start.size());
411 
412  // Work arrays
413  labelList hitMap(identity(start.size()));
414  pointField p0(start);
415  pointField p1(end);
416  List<pointIndexHit> intersectInfo(start.size());
417 
418  forAll(surfacesToTest, testI)
419  {
420  // Do synchronised call to all surfaces.
421  allSurfaces[surfacesToTest[testI]].findLineAny(p0, p1, intersectInfo);
422 
423  // Copy all hits into arguments, continue with misses
424  label newI = 0;
425  forAll(intersectInfo, i)
426  {
427  if (intersectInfo[i].hit())
428  {
429  hitInfo[hitMap[i]] = intersectInfo[i];
430  hitSurfaces[hitMap[i]] = testI;
431  }
432  else
433  {
434  if (i != newI)
435  {
436  hitMap[newI] = hitMap[i];
437  p0[newI] = p0[i];
438  p1[newI] = p1[i];
439  }
440  newI++;
441  }
442  }
443 
444  // All done? Note that this decision should be synchronised
445  if (newI == 0)
446  {
447  break;
448  }
449 
450  // Trim and continue
451  hitMap.setSize(newI);
452  p0.setSize(newI);
453  p1.setSize(newI);
454  intersectInfo.setSize(newI);
455  }
456 }
457 
458 
460 (
461  const PtrList<searchableSurface>& allSurfaces,
462  const labelList& surfacesToTest,
463  const pointField& start,
464  const pointField& end,
465  labelListList& hitSurfaces,
466  List<List<pointIndexHit> >& hitInfo
467 )
468 {
469  // Note: maybe move the single-surface all intersections test into
470  // searchable surface? Some of the tolerance issues might be
471  // lessened.
472 
473  // 2. Currently calling searchableSurface::findLine with start==end
474  // is expected to find no intersection. Problem if it does.
475 
476  hitSurfaces.setSize(start.size());
477  hitInfo.setSize(start.size());
478 
479  if (surfacesToTest.empty())
480  {
481  return;
482  }
483 
484  // Test first surface
485  allSurfaces[surfacesToTest[0]].findLineAll(start, end, hitInfo);
486 
487  // Set hitSurfaces and distance
488  List<scalarList> hitDistSqr(hitInfo.size());
489  forAll(hitInfo, pointI)
490  {
491  const List<pointIndexHit>& pHits = hitInfo[pointI];
492 
493  labelList& pSurfaces = hitSurfaces[pointI];
494  pSurfaces.setSize(pHits.size());
495  pSurfaces = 0;
496 
497  scalarList& pDistSqr = hitDistSqr[pointI];
498  pDistSqr.setSize(pHits.size());
499  forAll(pHits, i)
500  {
501  pDistSqr[i] = magSqr(pHits[i].hitPoint() - start[pointI]);
502  }
503  }
504 
505 
506  if (surfacesToTest.size() > 1)
507  {
508  // Small vector to increment start vector by to find next intersection
509  // along line. Constant factor added to make sure that start==end still
510  // ends iteration in findAllIntersections. Also SMALL is just slightly
511  // too small.
512  const vectorField smallVec
513  (
514  1E2*SMALL*(end-start)
515  + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
516  );
517 
518  // Tolerance used to check whether points are equal. Note: used to
519  // compare distance^2. Note that we use the maximum possible tolerance
520  // (reached at intersections close to the end point)
521  const scalarField mergeDist(2*mag(smallVec)*mag(end-start));
522 
523  // Test the other surfaces and merge (according to distance from start).
524  for (label testI = 1; testI < surfacesToTest.size(); testI++)
525  {
526  List<List<pointIndexHit> > surfHits;
527  allSurfaces[surfacesToTest[testI]].findLineAll
528  (
529  start,
530  end,
531  surfHits
532  );
533 
534  forAll(surfHits, pointI)
535  {
536  mergeHits
537  (
538  start[pointI], // Current segment
539  mergeDist[pointI],
540 
541  testI, // Surface and its hits
542  surfHits[pointI],
543 
544  hitSurfaces[pointI], // Merge into overall hit info
545  hitInfo[pointI],
546  hitDistSqr[pointI]
547  );
548  }
549  }
550  }
551 }
552 
553 
555 //void Foam::searchableSurfacesQueries::findNearestIntersection
556 //(
557 // const PtrList<searchableSurface>& allSurfaces,
558 // const labelList& surfacesToTest,
559 // const pointField& start,
560 // const pointField& end,
561 //
562 // labelList& surface1,
563 // List<pointIndexHit>& hit1,
564 // labelList& surface2,
565 // List<pointIndexHit>& hit2
566 //)
567 //{
568 // // 1. intersection from start to end
569 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
570 //
571 // // Initialize arguments
572 // surface1.setSize(start.size());
573 // surface1 = -1;
574 // hit1.setSize(start.size());
575 //
576 // // Current end of segment to test.
577 // pointField nearest(end);
578 // // Work array
579 // List<pointIndexHit> nearestInfo(start.size());
580 //
581 // forAll(surfacesToTest, testI)
582 // {
583 // // See if any intersection between start and current nearest
584 // allSurfaces[surfacesToTest[testI]].findLine
585 // (
586 // start,
587 // nearest,
588 // nearestInfo
589 // );
590 //
591 // forAll(nearestInfo, pointI)
592 // {
593 // if (nearestInfo[pointI].hit())
594 // {
595 // hit1[pointI] = nearestInfo[pointI];
596 // surface1[pointI] = testI;
597 // nearest[pointI] = hit1[pointI].hitPoint();
598 // }
599 // }
600 // }
601 //
602 //
603 // // 2. intersection from end to last intersection
604 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
605 //
606 // // Find the nearest intersection from end to start. Note that we
607 // // initialize to the first intersection (if any).
608 // surface2 = surface1;
609 // hit2 = hit1;
610 //
611 // // Set current end of segment to test.
612 // forAll(nearest, pointI)
613 // {
614 // if (hit1[pointI].hit())
615 // {
616 // nearest[pointI] = hit1[pointI].hitPoint();
617 // }
618 // else
619 // {
620 // // Disable testing by setting to end.
621 // nearest[pointI] = end[pointI];
622 // }
623 // }
624 //
625 // forAll(surfacesToTest, testI)
626 // {
627 // // See if any intersection between end and current nearest
628 // allSurfaces[surfacesToTest[i]].findLine(end, nearest, nearestInfo);
629 //
630 // forAll(nearestInfo, pointI)
631 // {
632 // if (nearestInfo[pointI].hit())
633 // {
634 // hit2[pointI] = nearestInfo[pointI];
635 // surface2[pointI] = testI;
636 // nearest[pointI] = hit2[pointI].hitPoint();
637 // }
638 // }
639 // }
640 //}
641 
642 
643 // Find nearest. Return -1 or nearest point
645 (
646  const PtrList<searchableSurface>& allSurfaces,
647  const labelList& surfacesToTest,
648  const pointField& samples,
649  const scalarField& nearestDistSqr,
650  labelList& nearestSurfaces,
651  List<pointIndexHit>& nearestInfo
652 )
653 {
654  // Initialise
655  nearestSurfaces.setSize(samples.size());
656  nearestSurfaces = -1;
657  nearestInfo.setSize(samples.size());
658 
659  // Work arrays
660  scalarField minDistSqr(nearestDistSqr);
661  List<pointIndexHit> hitInfo(samples.size());
662 
663  forAll(surfacesToTest, testI)
664  {
665  allSurfaces[surfacesToTest[testI]].findNearest
666  (
667  samples,
668  minDistSqr,
669  hitInfo
670  );
671 
672  // Update minDistSqr and arguments
673  forAll(hitInfo, pointI)
674  {
675  if (hitInfo[pointI].hit())
676  {
677  minDistSqr[pointI] = magSqr
678  (
679  hitInfo[pointI].hitPoint()
680  - samples[pointI]
681  );
682  nearestInfo[pointI] = hitInfo[pointI];
683  nearestSurfaces[pointI] = testI;
684  }
685  }
686  }
687 }
688 
689 
690 //- Calculate point which is on a set of surfaces.
692 (
693  const PtrList<searchableSurface>& allSurfaces,
694  const labelList& surfacesToTest,
695  const scalar initDistSqr,
696  const scalar convergenceDistSqr,
697  const point& start
698 )
699 {
700  // Get four starting points. Take these as the projection of the
701  // starting point onto the surfaces and the mid point
702  List<point> nearest(surfacesToTest.size()+1);
703 
704  point sumNearest = vector::zero;
705 
706  forAll(surfacesToTest, i)
707  {
708  pointIndexHit hit
709  (
710  tempFindNearest(allSurfaces[surfacesToTest[i]], start, initDistSqr)
711  );
712 
713  if (hit.hit())
714  {
715  nearest[i] = hit.hitPoint();
716  sumNearest += nearest[i];
717  }
718  else
719  {
721  (
722  "searchableSurfacesQueries::facesIntersection"
723  "(const labelList&, const scalar, const scalar, const point&)"
724  ) << "Did not find point within distance "
725  << initDistSqr << " of starting point " << start
726  << " on surface "
727  << allSurfaces[surfacesToTest[i]].IOobject::name()
728  << abort(FatalError);
729  }
730  }
731 
732  nearest[nearest.size()-1] = sumNearest / surfacesToTest.size();
733 
734 
735  // Get the sum of distances (initial evaluation)
736  List<scalar> nearestDist(nearest.size());
737 
738  forAll(nearestDist, i)
739  {
740  nearestDist[i] = sumDistSqr
741  (
742  allSurfaces,
743  surfacesToTest,
744  initDistSqr,
745  nearest[i]
746  );
747  }
748 
749 
750  //- Downhill Simplex method
751 
752  bool converged = morphTet
753  (
754  allSurfaces,
755  surfacesToTest,
756  initDistSqr,
757  convergenceDistSqr,
758  2000,
759  nearest,
760  nearestDist
761  );
762 
763 
765 
766  if (converged)
767  {
768  // Project nearest onto 0th surface.
769  intersection = tempFindNearest
770  (
771  allSurfaces[surfacesToTest[0]],
772  nearest[0],
773  nearestDist[0]
774  );
775  }
776 
777  //if (!intersection.hit())
778  //{
779  // // Restart
780  // scalar smallDist = Foam::sqr(convergenceDistSqr);
781  // nearest[0] = intersection.hitPoint();
782  // nearest[1] = nearest[0];
783  // nearest[1].x() += smallDist;
784  // nearest[2] = nearest[0];
785  // nearest[2].y() += smallDist;
786  // nearest[3] = nearest[0];
787  // nearest[3].z() += smallDist;
788  //
789  // forAll(nearestDist, i)
790  // {
791  // nearestDist[i] = sumDistSqr
792  // (
793  // surfacesToTest,
794  // initDistSqr,
795  // nearest[i]
796  // );
797  // }
798  //
799  // intersection = morphTet
800  // (
801  // allSurfaces,
802  // surfacesToTest,
803  // initDistSqr,
804  // convergenceDistSqr,
805  // 1000,
806  // nearest,
807  // nearestDist
808  // );
809  //}
810 
811  return intersection;
812 }
813 
814 
815 
816 // ************************ vim: set sw=4 sts=4 et: ************************ //