FreeFOAM The Cross-Platform CFD Toolkit
searchableSurfaceCollection.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 
28 #include <OpenFOAM/SortableList.H>
29 #include <OpenFOAM/Time.H>
30 #include <OpenFOAM/ListOps.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(searchableSurfaceCollection, 0);
39 (
40  searchableSurface,
41  searchableSurfaceCollection,
42  dict
43 );
44 
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void Foam::searchableSurfaceCollection::findNearest
51 (
52  const pointField& samples,
53  scalarField& minDistSqr,
54  List<pointIndexHit>& nearestInfo,
55  labelList& nearestSurf
56 ) const
57 {
58  // Initialise
59  nearestInfo.setSize(samples.size());
60  nearestInfo = pointIndexHit();
61  nearestSurf.setSize(samples.size());
62  nearestSurf = -1;
63 
64  List<pointIndexHit> hitInfo(samples.size());
65 
66  const scalarField localMinDistSqr(samples.size(), GREAT);
67 
68  forAll(subGeom_, surfI)
69  {
70  subGeom_[surfI].findNearest
71  (
72  cmptDivide // Transform then divide
73  (
74  transform_[surfI].localPosition(samples),
75  scale_[surfI]
76  ),
77  localMinDistSqr,
78  hitInfo
79  );
80 
81  forAll(hitInfo, pointI)
82  {
83  if (hitInfo[pointI].hit())
84  {
85  // Rework back into global coordinate sys. Multiply then
86  // transform
87  point globalPt = transform_[surfI].globalPosition
88  (
90  (
91  hitInfo[pointI].rawPoint(),
92  scale_[surfI]
93  )
94  );
95 
96  scalar distSqr = magSqr(globalPt - samples[pointI]);
97 
98  if (distSqr < minDistSqr[pointI])
99  {
100  minDistSqr[pointI] = distSqr;
101  nearestInfo[pointI].setPoint(globalPt);
102  nearestInfo[pointI].setHit();
103  nearestInfo[pointI].setIndex
104  (
105  hitInfo[pointI].index()
106  + indexOffset_[surfI]
107  );
108  nearestSurf[pointI] = surfI;
109  }
110  }
111  }
112  }
113 }
114 
115 
116 // Sort hits into per-surface bins. Misses are rejected. Maintains map back
117 // to position
118 void Foam::searchableSurfaceCollection::sortHits
119 (
120  const List<pointIndexHit>& info,
121  List<List<pointIndexHit> >& surfInfo,
122  labelListList& infoMap
123 ) const
124 {
125  // Count hits per surface.
126  labelList nHits(subGeom_.size(), 0);
127 
128  forAll(info, pointI)
129  {
130  if (info[pointI].hit())
131  {
132  label index = info[pointI].index();
133  label surfI = findLower(indexOffset_, index+1);
134  nHits[surfI]++;
135  }
136  }
137 
138  // Per surface the hit
139  surfInfo.setSize(subGeom_.size());
140  // Per surface the original position
141  infoMap.setSize(subGeom_.size());
142 
143  forAll(surfInfo, surfI)
144  {
145  surfInfo[surfI].setSize(nHits[surfI]);
146  infoMap[surfI].setSize(nHits[surfI]);
147  }
148  nHits = 0;
149 
150  forAll(info, pointI)
151  {
152  if (info[pointI].hit())
153  {
154  label index = info[pointI].index();
155  label surfI = findLower(indexOffset_, index+1);
156 
157  // Store for correct surface and adapt indices back to local
158  // ones
159  label localI = nHits[surfI]++;
160  surfInfo[surfI][localI] = pointIndexHit
161  (
162  info[pointI].hit(),
163  info[pointI].rawPoint(),
164  index-indexOffset_[surfI]
165  );
166  infoMap[surfI][localI] = pointI;
167  }
168  }
169 }
170 
171 
172 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 
174 Foam::searchableSurfaceCollection::searchableSurfaceCollection
175 (
176  const IOobject& io,
177  const dictionary& dict
178 )
179 :
180  searchableSurface(io),
181  instance_(dict.size()),
182  scale_(dict.size()),
183  transform_(dict.size()),
184  subGeom_(dict.size()),
185  mergeSubRegions_(dict.lookup("mergeSubRegions")),
186  indexOffset_(dict.size()+1)
187 {
188  Info<< "SearchableCollection : " << name() << endl;
189 
190  label surfI = 0;
191  label startIndex = 0;
192  forAllConstIter(dictionary, dict, iter)
193  {
194  if (dict.isDict(iter().keyword()))
195  {
196  instance_[surfI] = iter().keyword();
197 
198  const dictionary& subDict = dict.subDict(instance_[surfI]);
199 
200  scale_[surfI] = subDict.lookup("scale");
201  transform_.set
202  (
203  surfI,
205  (
206  "",
207  subDict.subDict("transform")
208  )
209  );
210 
211  const word subGeomName(subDict.lookup("surface"));
212  //Pout<< "Trying to find " << subGeomName << endl;
213 
214  const searchableSurface& s =
215  io.db().lookupObject<searchableSurface>(subGeomName);
216 
217  // I don't know yet how to handle the globalSize combined with
218  // regionOffset. Would cause non-consecutive indices locally
219  // if all indices offset by globalSize() of the local region...
220  if (s.size() != s.globalSize())
221  {
223  (
224  "searchableSurfaceCollection::searchableSurfaceCollection"
225  "(const IOobject&, const dictionary&)"
226  ) << "Cannot use a distributed surface in a collection."
227  << exit(FatalError);
228  }
229 
230  subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
231 
232  indexOffset_[surfI] = startIndex;
233  startIndex += subGeom_[surfI].size();
234 
235  Info<< " instance : " << instance_[surfI] << endl;
236  Info<< " surface : " << s.name() << endl;
237  Info<< " scale : " << scale_[surfI] << endl;
238  Info<< " coordsys : " << transform_[surfI] << endl;
239 
240  surfI++;
241  }
242  }
243  indexOffset_[surfI] = startIndex;
244 
245  instance_.setSize(surfI);
246  scale_.setSize(surfI);
247  transform_.setSize(surfI);
248  subGeom_.setSize(surfI);
249  indexOffset_.setSize(surfI+1);
250 }
251 
252 
253 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
254 
256 {}
257 
258 
259 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
260 
262 {
263  if (regions_.size() == 0)
264  {
265  regionOffset_.setSize(subGeom_.size());
266 
267  DynamicList<word> allRegions;
268  forAll(subGeom_, surfI)
269  {
270  regionOffset_[surfI] = allRegions.size();
271 
272  if (mergeSubRegions_)
273  {
274  // Single name regardless how many regions subsurface has
275  allRegions.append(instance_[surfI] + "_" + Foam::name(surfI));
276  }
277  else
278  {
279  const wordList& subRegions = subGeom_[surfI].regions();
280 
281  forAll(subRegions, i)
282  {
283  allRegions.append(instance_[surfI] + "_" + subRegions[i]);
284  }
285  }
286  }
287  regions_.transfer(allRegions.shrink());
288  }
289  return regions_;
290 }
291 
292 
294 {
295  return indexOffset_[indexOffset_.size()-1];
296 }
297 
298 
300 {
301  // Get overall size
302  pointField coords(size());
303 
304  // Append individual coordinates
305  label coordI = 0;
306 
307  forAll(subGeom_, surfI)
308  {
309  const pointField subCoords = subGeom_[surfI].coordinates();
310 
311  forAll(subCoords, i)
312  {
313  coords[coordI++] = transform_[surfI].globalPosition
314  (
316  (
317  subCoords[i],
318  scale_[surfI]
319  )
320  );
321  }
322  }
323 
324  return coords;
325 }
326 
327 
328 void Foam::searchableSurfaceCollection::findNearest
329 (
330  const pointField& samples,
331  const scalarField& nearestDistSqr,
332  List<pointIndexHit>& nearestInfo
333 ) const
334 {
335  // How to scale distance?
336  scalarField minDistSqr(nearestDistSqr);
337 
338  labelList nearestSurf;
339  findNearest
340  (
341  samples,
342  minDistSqr,
343  nearestInfo,
344  nearestSurf
345  );
346 }
347 
348 
350 (
351  const pointField& start,
352  const pointField& end,
353  List<pointIndexHit>& info
354 ) const
355 {
356  info.setSize(start.size());
357  info = pointIndexHit();
358 
359  // Current nearest (to start) intersection
360  pointField nearest(end);
361 
362  List<pointIndexHit> hitInfo(start.size());
363 
364  forAll(subGeom_, surfI)
365  {
366  // Starting point
368  (
369  transform_[surfI].localPosition
370  (
371  start
372  ),
373  scale_[surfI]
374  );
375 
376  // Current best end point
378  (
379  transform_[surfI].localPosition
380  (
381  nearest
382  ),
383  scale_[surfI]
384  );
385 
386  subGeom_[surfI].findLine(e0, e1, hitInfo);
387 
388  forAll(hitInfo, pointI)
389  {
390  if (hitInfo[pointI].hit())
391  {
392  // Transform back to global coordinate sys.
393  nearest[pointI] = transform_[surfI].globalPosition
394  (
396  (
397  hitInfo[pointI].rawPoint(),
398  scale_[surfI]
399  )
400  );
401  info[pointI] = hitInfo[pointI];
402  info[pointI].rawPoint() = nearest[pointI];
403  info[pointI].setIndex
404  (
405  hitInfo[pointI].index()
406  + indexOffset_[surfI]
407  );
408  }
409  }
410  }
411 
412 
413  // Debug check
414  if (false)
415  {
416  forAll(info, pointI)
417  {
418  if (info[pointI].hit())
419  {
420  vector n(end[pointI] - start[pointI]);
421  scalar magN = mag(n);
422 
423  if (magN > SMALL)
424  {
425  n /= mag(n);
426 
427  scalar s = ((info[pointI].rawPoint()-start[pointI])&n);
428 
429  if (s < 0 || s > 1)
430  {
432  (
433  "searchableSurfaceCollection::findLine(..)"
434  ) << "point:" << info[pointI]
435  << " s:" << s
436  << " outside vector "
437  << " start:" << start[pointI]
438  << " end:" << end[pointI]
439  << abort(FatalError);
440  }
441  }
442  }
443  }
444  }
445 }
446 
447 
449 (
450  const pointField& start,
451  const pointField& end,
452  List<pointIndexHit>& info
453 ) const
454 {
455  // To be done ...
456  findLine(start, end, info);
457 }
458 
459 
461 (
462  const pointField& start,
463  const pointField& end,
464  List<List<pointIndexHit> >& info
465 ) const
466 {
467  // To be done. Assume for now only one intersection.
468  List<pointIndexHit> nearestInfo;
469  findLine(start, end, nearestInfo);
470 
471  info.setSize(start.size());
472  forAll(info, pointI)
473  {
474  if (nearestInfo[pointI].hit())
475  {
476  info[pointI].setSize(1);
477  info[pointI][0] = nearestInfo[pointI];
478  }
479  else
480  {
481  info[pointI].clear();
482  }
483  }
484 }
485 
486 
488 (
489  const List<pointIndexHit>& info,
490  labelList& region
491 ) const
492 {
493  if (subGeom_.size() == 0)
494  {}
495  else if (subGeom_.size() == 1)
496  {
497  if (mergeSubRegions_)
498  {
499  region.setSize(info.size());
500  region = regionOffset_[0];
501  }
502  else
503  {
504  subGeom_[0].getRegion(info, region);
505  }
506  }
507  else
508  {
509  // Multiple surfaces. Sort by surface.
510 
511  // Per surface the hit
512  List<List<pointIndexHit> > surfInfo;
513  // Per surface the original position
514  List<List<label> > infoMap;
515  sortHits(info, surfInfo, infoMap);
516 
517  region.setSize(info.size());
518  region = -1;
519 
520  // Do region tests
521 
522  if (mergeSubRegions_)
523  {
524  // Actually no need for surfInfo. Just take region for surface.
525  forAll(infoMap, surfI)
526  {
527  const labelList& map = infoMap[surfI];
528  forAll(map, i)
529  {
530  region[map[i]] = regionOffset_[surfI];
531  }
532  }
533  }
534  else
535  {
536  forAll(infoMap, surfI)
537  {
538  labelList surfRegion;
539  subGeom_[surfI].getRegion(surfInfo[surfI], surfRegion);
540 
541  const labelList& map = infoMap[surfI];
542  forAll(map, i)
543  {
544  region[map[i]] = regionOffset_[surfI] + surfRegion[i];
545  }
546  }
547  }
548  }
549 }
550 
551 
553 (
554  const List<pointIndexHit>& info,
556 ) const
557 {
558  if (subGeom_.size() == 0)
559  {}
560  else if (subGeom_.size() == 1)
561  {
562  subGeom_[0].getNormal(info, normal);
563  }
564  else
565  {
566  // Multiple surfaces. Sort by surface.
567 
568  // Per surface the hit
569  List<List<pointIndexHit> > surfInfo;
570  // Per surface the original position
571  List<List<label> > infoMap;
572  sortHits(info, surfInfo, infoMap);
573 
574  normal.setSize(info.size());
575 
576  // Do region tests
577  forAll(surfInfo, surfI)
578  {
579  vectorField surfNormal;
580  subGeom_[surfI].getNormal(surfInfo[surfI], surfNormal);
581 
582  const labelList& map = infoMap[surfI];
583  forAll(map, i)
584  {
585  normal[map[i]] = surfNormal[i];
586  }
587  }
588  }
589 }
590 
591 
593 (
594  const pointField& points,
595  List<volumeType>& volType
596 ) const
597 {
599  (
600  "searchableSurfaceCollection::getVolumeType(const pointField&"
601  ", List<volumeType>&) const"
602  ) << "Volume type not supported for collection."
603  << exit(FatalError);
604 }
605 
606 
608 (
609  const List<treeBoundBox>& bbs,
610  const bool keepNonLocal,
611  autoPtr<mapDistribute>& faceMap,
612  autoPtr<mapDistribute>& pointMap
613 )
614 {
615  forAll(subGeom_, surfI)
616  {
617  // Note:Tranform the bounding boxes? Something like
618  // pointField bbPoints =
619  // cmptDivide
620  // (
621  // transform_[surfI].localPosition
622  // (
623  // bbs[i].points()
624  // ),
625  // scale_[surfI]
626  // );
627  // treeBoundBox newBb(bbPoints);
628 
629  // Note: what to do with faceMap, pointMap from multiple surfaces?
630  subGeom_[surfI].distribute
631  (
632  bbs,
633  keepNonLocal,
634  faceMap,
635  pointMap
636  );
637  }
638 }
639 
640 
642 {
643  forAll(subGeom_, surfI)
644  {
645  subGeom_[surfI].setField
646  (
647  static_cast<const labelList&>
648  (
650  (
651  values,
652  subGeom_[surfI].size(),
653  indexOffset_[surfI]
654  )
655  )
656  );
657  }
658 }
659 
660 
662 (
663  const List<pointIndexHit>& info,
664  labelList& values
665 ) const
666 {
667  if (subGeom_.size() == 0)
668  {}
669  else if (subGeom_.size() == 1)
670  {
671  subGeom_[0].getField(info, values);
672  }
673  else
674  {
675  // Multiple surfaces. Sort by surface.
676 
677  // Per surface the hit
678  List<List<pointIndexHit> > surfInfo;
679  // Per surface the original position
680  List<List<label> > infoMap;
681  sortHits(info, surfInfo, infoMap);
682 
683  // Do surface tests
684  forAll(surfInfo, surfI)
685  {
686  labelList surfValues;
687  subGeom_[surfI].getField(surfInfo[surfI], surfValues);
688 
689  if (surfValues.size())
690  {
691  // Size values only when we have a surface that supports it.
692  values.setSize(info.size());
693 
694  const labelList& map = infoMap[surfI];
695  forAll(map, i)
696  {
697  values[map[i]] = surfValues[i];
698  }
699  }
700  }
701  }
702 }
703 
704 
705 // ************************ vim: set sw=4 sts=4 et: ************************ //