distributedTriSurfaceMesh.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-2017 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 \*---------------------------------------------------------------------------*/
25 
27 #include "mapDistribute.H"
28 #include "Random.H"
30 #include "triangleFuncs.H"
31 #include "matchPoints.H"
32 #include "globalIndex.H"
33 #include "Time.H"
34 
35 #include "IFstream.H"
36 #include "decompositionMethod.H"
37 #include "geomDecomp.H"
38 #include "vectorList.H"
39 #include "PackedBoolList.H"
40 #include "PatchTools.H"
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 namespace Foam
45 {
46  defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
48  (
49  searchableSurface,
50  distributedTriSurfaceMesh,
51  dict
52  );
53 
54  template<>
55  const char* Foam::NamedEnum
56  <
58  3
59  >::names[] =
60  {
61  "follow",
62  "independent",
63  "frozen"
64  };
65 }
66 
67 
70 
71 
72 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
73 
74 // Read my additional data from the dictionary
75 bool Foam::distributedTriSurfaceMesh::read()
76 {
77  // Get bb of all domains.
78  procBb_.setSize(Pstream::nProcs());
79 
80  procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
81  Pstream::gatherList(procBb_);
82  Pstream::scatterList(procBb_);
83 
84  // Distribution type
85  distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
86 
87  // Merge distance
88  mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
89 
90  return true;
91 }
92 
93 
94 // Is segment fully local?
95 bool Foam::distributedTriSurfaceMesh::isLocal
96 (
97  const List<treeBoundBox>& myBbs,
98  const point& start,
99  const point& end
100 )
101 {
102  forAll(myBbs, bbI)
103  {
104  if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
105  {
106  return true;
107  }
108  }
109  return false;
110 }
111 
112 
113 //void Foam::distributedTriSurfaceMesh::splitSegment
114 //(
115 // const label segmentI,
116 // const point& start,
117 // const point& end,
118 // const treeBoundBox& bb,
119 //
120 // DynamicList<segment>& allSegments,
121 // DynamicList<label>& allSegmentMap,
122 // DynamicList<label> sendMap
123 //) const
124 //{
125 // // Work points
126 // point clipPt0, clipPt1;
127 //
128 // if (bb.contains(start))
129 // {
130 // // start within, trim end to bb
131 // bool clipped = bb.intersects(end, start, clipPt0);
132 //
133 // if (clipped)
134 // {
135 // // segment from start to clippedStart passes
136 // // through proc.
137 // sendMap[proci].append(allSegments.size());
138 // allSegmentMap.append(segmentI);
139 // allSegments.append(segment(start, clipPt0));
140 // }
141 // }
142 // else if (bb.contains(end))
143 // {
144 // // end within, trim start to bb
145 // bool clipped = bb.intersects(start, end, clipPt0);
146 //
147 // if (clipped)
148 // {
149 // sendMap[proci].append(allSegments.size());
150 // allSegmentMap.append(segmentI);
151 // allSegments.append(segment(clipPt0, end));
152 // }
153 // }
154 // else
155 // {
156 // // trim both
157 // bool clippedStart = bb.intersects(start, end, clipPt0);
158 //
159 // if (clippedStart)
160 // {
161 // bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
162 //
163 // if (clippedEnd)
164 // {
165 // // middle part of segment passes through proc.
166 // sendMap[proci].append(allSegments.size());
167 // allSegmentMap.append(segmentI);
168 // allSegments.append(segment(clipPt0, clipPt1));
169 // }
170 // }
171 // }
172 //}
173 
174 
175 void Foam::distributedTriSurfaceMesh::distributeSegment
176 (
177  const label segmentI,
178  const point& start,
179  const point& end,
180 
181  DynamicList<segment>& allSegments,
182  DynamicList<label>& allSegmentMap,
183  List<DynamicList<label>>& sendMap
184 ) const
185 {
186  // 1. Fully local already handled outside. Note: retest is cheap.
187  if (isLocal(procBb_[Pstream::myProcNo()], start, end))
188  {
189  return;
190  }
191 
192 
193  // 2. If fully inside one other processor, then only need to send
194  // to that one processor even if it intersects another. Rare occurrence
195  // but cheap to test.
196  forAll(procBb_, proci)
197  {
198  if (proci != Pstream::myProcNo())
199  {
200  const List<treeBoundBox>& bbs = procBb_[proci];
201 
202  if (isLocal(bbs, start, end))
203  {
204  sendMap[proci].append(allSegments.size());
205  allSegmentMap.append(segmentI);
206  allSegments.append(segment(start, end));
207  return;
208  }
209  }
210  }
211 
212 
213  // 3. If not contained in single processor send to all intersecting
214  // processors.
215  forAll(procBb_, proci)
216  {
217  const List<treeBoundBox>& bbs = procBb_[proci];
218 
219  forAll(bbs, bbI)
220  {
221  const treeBoundBox& bb = bbs[bbI];
222 
223  // Scheme a: any processor that intersects the segment gets
224  // the segment.
225 
226  // Intersection point
227  point clipPt;
228 
229  if (bb.intersects(start, end, clipPt))
230  {
231  sendMap[proci].append(allSegments.size());
232  allSegmentMap.append(segmentI);
233  allSegments.append(segment(start, end));
234  }
235 
236  // Alternative: any processor only gets clipped bit of
237  // segment. This gives small problems with additional
238  // truncation errors.
239  //splitSegment
240  //(
241  // segmentI,
242  // start,
243  // end,
244  // bb,
245  //
246  // allSegments,
247  // allSegmentMap,
248  // sendMap[proci]
249  //);
250  }
251  }
252 }
253 
254 
256 Foam::distributedTriSurfaceMesh::distributeSegments
257 (
258  const pointField& start,
259  const pointField& end,
260 
261  List<segment>& allSegments,
262  labelList& allSegmentMap
263 ) const
264 {
265  // Determine send map
266  // ~~~~~~~~~~~~~~~~~~
267 
268  labelListList sendMap(Pstream::nProcs());
269 
270  {
271  // Since intersection test is quite expensive compared to memory
272  // allocation we use DynamicList to immediately store the segment
273  // in the correct bin.
274 
275  // Segments to test
276  DynamicList<segment> dynAllSegments(start.size());
277  // Original index of segment
278  DynamicList<label> dynAllSegmentMap(start.size());
279  // Per processor indices into allSegments to send
281 
282  forAll(start, segmentI)
283  {
284  distributeSegment
285  (
286  segmentI,
287  start[segmentI],
288  end[segmentI],
289 
290  dynAllSegments,
291  dynAllSegmentMap,
292  dynSendMap
293  );
294  }
295 
296  // Convert dynamicList to labelList
297  sendMap.setSize(Pstream::nProcs());
298  forAll(sendMap, proci)
299  {
300  dynSendMap[proci].shrink();
301  sendMap[proci].transfer(dynSendMap[proci]);
302  }
303 
304  allSegments.transfer(dynAllSegments.shrink());
305  allSegmentMap.transfer(dynAllSegmentMap.shrink());
306  }
307 
308 
309  // Send over how many I need to receive.
310  labelListList sendSizes(Pstream::nProcs());
311  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
312  forAll(sendMap, proci)
313  {
314  sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
315  }
316  Pstream::gatherList(sendSizes);
317  Pstream::scatterList(sendSizes);
318 
319 
320  // Determine order of receiving
321  labelListList constructMap(Pstream::nProcs());
322 
323  // My local segments first
324  constructMap[Pstream::myProcNo()] = identity
325  (
326  sendMap[Pstream::myProcNo()].size()
327  );
328 
329  label segmentI = constructMap[Pstream::myProcNo()].size();
330  forAll(constructMap, proci)
331  {
332  if (proci != Pstream::myProcNo())
333  {
334  // What I need to receive is what other processor is sending to me.
335  label nRecv = sendSizes[proci][Pstream::myProcNo()];
336  constructMap[proci].setSize(nRecv);
337 
338  for (label i = 0; i < nRecv; i++)
339  {
340  constructMap[proci][i] = segmentI++;
341  }
342  }
343  }
344 
346  (
347  new mapDistribute
348  (
349  segmentI, // size after construction
350  sendMap.xfer(),
351  constructMap.xfer()
352  )
353  );
354 }
355 
356 
357 void Foam::distributedTriSurfaceMesh::findLine
358 (
359  const bool nearestIntersection,
360  const pointField& start,
361  const pointField& end,
362  List<pointIndexHit>& info
363 ) const
364 {
365  const indexedOctree<treeDataTriSurface>& octree = tree();
366 
367  // Initialise
368  info.setSize(start.size());
369  forAll(info, i)
370  {
371  info[i].setMiss();
372  }
373 
374  if (!Pstream::parRun())
375  {
376  forAll(start, i)
377  {
378  if (nearestIntersection)
379  {
380  info[i] = octree.findLine(start[i], end[i]);
381  }
382  else
383  {
384  info[i] = octree.findLineAny(start[i], end[i]);
385  }
386  }
387  }
388  else
389  {
390  // Important:force synchronised construction of indexing
391  const globalIndex& triIndexer = globalTris();
392 
393 
394  // Do any local queries
395  // ~~~~~~~~~~~~~~~~~~~~
396 
397  label nLocal = 0;
398 
399  forAll(start, i)
400  {
401  if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
402  {
403  if (nearestIntersection)
404  {
405  info[i] = octree.findLine(start[i], end[i]);
406  }
407  else
408  {
409  info[i] = octree.findLineAny(start[i], end[i]);
410  }
411 
412  if (info[i].hit())
413  {
414  info[i].setIndex(triIndexer.toGlobal(info[i].index()));
415  }
416  nLocal++;
417  }
418  }
419 
420 
421  if
422  (
423  returnReduce(nLocal, sumOp<label>())
424  < returnReduce(start.size(), sumOp<label>())
425  )
426  {
427  // Not all can be resolved locally. Build segments and map,
428  // send over segments, do intersections, send back and merge.
429 
430 
431  // Construct queries (segments)
432  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433 
434  // Segments to test
435  List<segment> allSegments(start.size());
436  // Original index of segment
437  labelList allSegmentMap(start.size());
438 
439  const autoPtr<mapDistribute> mapPtr
440  (
441  distributeSegments
442  (
443  start,
444  end,
445  allSegments,
446  allSegmentMap
447  )
448  );
449  const mapDistribute& map = mapPtr();
450 
451  label nOldAllSegments = allSegments.size();
452 
453 
454  // Exchange the segments
455  // ~~~~~~~~~~~~~~~~~~~~~
456 
457  map.distribute(allSegments);
458 
459 
460  // Do tests I need to do
461  // ~~~~~~~~~~~~~~~~~~~~~
462 
463  // Intersections
464  List<pointIndexHit> intersections(allSegments.size());
465 
466  forAll(allSegments, i)
467  {
468  if (nearestIntersection)
469  {
470  intersections[i] = octree.findLine
471  (
472  allSegments[i].first(),
473  allSegments[i].second()
474  );
475  }
476  else
477  {
478  intersections[i] = octree.findLineAny
479  (
480  allSegments[i].first(),
481  allSegments[i].second()
482  );
483  }
484 
485  // Convert triangle index to global numbering
486  if (intersections[i].hit())
487  {
488  intersections[i].setIndex
489  (
490  triIndexer.toGlobal(intersections[i].index())
491  );
492  }
493  }
494 
495 
496  // Exchange the intersections (opposite to segments)
497  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
498 
499  map.reverseDistribute(nOldAllSegments, intersections);
500 
501 
502  // Extract the hits
503  // ~~~~~~~~~~~~~~~~
504 
505  forAll(intersections, i)
506  {
507  const pointIndexHit& allInfo = intersections[i];
508  label segmentI = allSegmentMap[i];
509  pointIndexHit& hitInfo = info[segmentI];
510 
511  if (allInfo.hit())
512  {
513  if (!hitInfo.hit())
514  {
515  // No intersection yet so take this one
516  hitInfo = allInfo;
517  }
518  else if (nearestIntersection)
519  {
520  // Nearest intersection
521  if
522  (
523  magSqr(allInfo.hitPoint()-start[segmentI])
524  < magSqr(hitInfo.hitPoint()-start[segmentI])
525  )
526  {
527  hitInfo = allInfo;
528  }
529  }
530  }
531  }
532  }
533  }
534 }
535 
536 
537 // Exchanges indices to the processor they come from.
538 // - calculates exchange map
539 // - uses map to calculate local triangle index
541 Foam::distributedTriSurfaceMesh::calcLocalQueries
542 (
543  const List<pointIndexHit>& info,
544  labelList& triangleIndex
545 ) const
546 {
547  triangleIndex.setSize(info.size());
548 
549  const globalIndex& triIndexer = globalTris();
550 
551 
552  // Determine send map
553  // ~~~~~~~~~~~~~~~~~~
554 
555  // Since determining which processor the query should go to is
556  // cheap we do a multi-pass algorithm to save some memory temporarily.
557 
558  // 1. Count
559  labelList nSend(Pstream::nProcs(), 0);
560 
561  forAll(info, i)
562  {
563  if (info[i].hit())
564  {
565  label proci = triIndexer.whichProcID(info[i].index());
566  nSend[proci]++;
567  }
568  }
569 
570  // 2. Size sendMap
571  labelListList sendMap(Pstream::nProcs());
572  forAll(nSend, proci)
573  {
574  sendMap[proci].setSize(nSend[proci]);
575  nSend[proci] = 0;
576  }
577 
578  // 3. Fill sendMap
579  forAll(info, i)
580  {
581  if (info[i].hit())
582  {
583  label proci = triIndexer.whichProcID(info[i].index());
584  triangleIndex[i] = triIndexer.toLocal(proci, info[i].index());
585  sendMap[proci][nSend[proci]++] = i;
586  }
587  else
588  {
589  triangleIndex[i] = -1;
590  }
591  }
592 
593 
594  // Send over how many I need to receive
595  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
596 
597  labelListList sendSizes(Pstream::nProcs());
598  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
599  forAll(sendMap, proci)
600  {
601  sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
602  }
603  Pstream::gatherList(sendSizes);
604  Pstream::scatterList(sendSizes);
605 
606 
607  // Determine receive map
608  // ~~~~~~~~~~~~~~~~~~~~~
609 
610  labelListList constructMap(Pstream::nProcs());
611 
612  // My local segments first
613  constructMap[Pstream::myProcNo()] = identity
614  (
615  sendMap[Pstream::myProcNo()].size()
616  );
617 
618  label segmentI = constructMap[Pstream::myProcNo()].size();
619  forAll(constructMap, proci)
620  {
621  if (proci != Pstream::myProcNo())
622  {
623  // What I need to receive is what other processor is sending to me.
624  label nRecv = sendSizes[proci][Pstream::myProcNo()];
625  constructMap[proci].setSize(nRecv);
626 
627  for (label i = 0; i < nRecv; i++)
628  {
629  constructMap[proci][i] = segmentI++;
630  }
631  }
632  }
633 
634 
635  // Pack into distribution map
636  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
637 
639  (
640  new mapDistribute
641  (
642  segmentI, // size after construction
643  sendMap.xfer(),
644  constructMap.xfer()
645  )
646  );
647  const mapDistribute& map = mapPtr();
648 
649 
650  // Send over queries
651  // ~~~~~~~~~~~~~~~~~
652 
653  map.distribute(triangleIndex);
654 
655 
656  return mapPtr;
657 }
658 
659 
660 Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
661 (
662  const point& centre,
663  const scalar radiusSqr,
664  boolList& overlaps
665 ) const
666 {
667  overlaps = false;
668  label nOverlaps = 0;
669 
670  forAll(procBb_, proci)
671  {
672  const List<treeBoundBox>& bbs = procBb_[proci];
673 
674  forAll(bbs, bbI)
675  {
676  if (bbs[bbI].overlaps(centre, radiusSqr))
677  {
678  overlaps[proci] = true;
679  nOverlaps++;
680  break;
681  }
682  }
683  }
684  return nOverlaps;
685 }
686 
687 
688 // Generate queries for parallel distance calculation
689 // - calculates exchange map
690 // - uses map to exchange points and radius
692 Foam::distributedTriSurfaceMesh::calcLocalQueries
693 (
694  const pointField& centres,
695  const scalarField& radiusSqr,
696 
697  pointField& allCentres,
698  scalarField& allRadiusSqr,
699  labelList& allSegmentMap
700 ) const
701 {
702  // Determine queries
703  // ~~~~~~~~~~~~~~~~~
704 
705  labelListList sendMap(Pstream::nProcs());
706 
707  {
708  // Queries
709  DynamicList<point> dynAllCentres(centres.size());
710  DynamicList<scalar> dynAllRadiusSqr(centres.size());
711  // Original index of segment
712  DynamicList<label> dynAllSegmentMap(centres.size());
713  // Per processor indices into allSegments to send
715 
716  // Work array - whether processor bb overlaps the bounding sphere.
717  boolList procBbOverlaps(Pstream::nProcs());
718 
719  forAll(centres, centreI)
720  {
721  // Find the processor this sample+radius overlaps.
722  calcOverlappingProcs
723  (
724  centres[centreI],
725  radiusSqr[centreI],
726  procBbOverlaps
727  );
728 
729  forAll(procBbOverlaps, proci)
730  {
731  if (proci != Pstream::myProcNo() && procBbOverlaps[proci])
732  {
733  dynSendMap[proci].append(dynAllCentres.size());
734  dynAllSegmentMap.append(centreI);
735  dynAllCentres.append(centres[centreI]);
736  dynAllRadiusSqr.append(radiusSqr[centreI]);
737  }
738  }
739  }
740 
741  // Convert dynamicList to labelList
742  sendMap.setSize(Pstream::nProcs());
743  forAll(sendMap, proci)
744  {
745  dynSendMap[proci].shrink();
746  sendMap[proci].transfer(dynSendMap[proci]);
747  }
748 
749  allCentres.transfer(dynAllCentres.shrink());
750  allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
751  allSegmentMap.transfer(dynAllSegmentMap.shrink());
752  }
753 
754 
755  // Send over how many I need to receive.
756  labelListList sendSizes(Pstream::nProcs());
757  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
758  forAll(sendMap, proci)
759  {
760  sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
761  }
762  Pstream::gatherList(sendSizes);
763  Pstream::scatterList(sendSizes);
764 
765 
766  // Determine order of receiving
767  labelListList constructMap(Pstream::nProcs());
768 
769  // My local segments first
770  constructMap[Pstream::myProcNo()] = identity
771  (
772  sendMap[Pstream::myProcNo()].size()
773  );
774 
775  label segmentI = constructMap[Pstream::myProcNo()].size();
776  forAll(constructMap, proci)
777  {
778  if (proci != Pstream::myProcNo())
779  {
780  // What I need to receive is what other processor is sending to me.
781  label nRecv = sendSizes[proci][Pstream::myProcNo()];
782  constructMap[proci].setSize(nRecv);
783 
784  for (label i = 0; i < nRecv; i++)
785  {
786  constructMap[proci][i] = segmentI++;
787  }
788  }
789  }
790 
792  (
793  new mapDistribute
794  (
795  segmentI, // size after construction
796  sendMap.xfer(),
797  constructMap.xfer()
798  )
799  );
800  return mapPtr;
801 }
802 
803 
804 // Find bounding boxes that guarantee a more or less uniform distribution
805 // of triangles. Decomposition in here is only used to get the bounding
806 // boxes, actual decomposition is done later on.
807 // Returns a per processor a list of bounding boxes that most accurately
808 // describe the shape. For now just a single bounding box per processor but
809 // optimisation might be to determine a better fitting shape.
811 Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
812 (
813  const triSurface& s
814 )
815 {
816  if (!decomposer_.valid())
817  {
818  // Use current decomposer.
819  // Note: or always use hierarchical?
820  IOdictionary decomposeDict
821  (
822  IOobject
823  (
824  "decomposeParDict",
829  false
830  )
831  );
832  decomposer_ = decompositionMethod::New(decomposeDict);
833 
834  if (!decomposer_().parallelAware())
835  {
837  << "The decomposition method " << decomposer_().typeName
838  << " does not decompose in parallel."
839  << " Please choose one that does." << exit(FatalError);
840  }
841 
842  if (!isA<geomDecomp>(decomposer_()))
843  {
845  << "The decomposition method " << decomposer_().typeName
846  << " is not a geometric decomposition method." << endl
847  << "Only geometric decomposition methods are currently"
848  << " supported."
849  << exit(FatalError);
850  }
851  }
852 
853  // Do decomposition according to triangle centre
854  pointField triCentres(s.size());
855  forAll(s, triI)
856  {
857  triCentres[triI] = s[triI].centre(s.points());
858  }
859 
860 
861  geomDecomp& decomposer = refCast<geomDecomp>(decomposer_());
862 
863  // Do the actual decomposition
864  labelList distribution(decomposer.decompose(triCentres));
865 
866  // Find bounding box for all triangles on new distribution.
867 
868  // Initialise to inverted box (VGREAT, -VGREAT)
870  forAll(bbs, proci)
871  {
872  bbs[proci].setSize(1);
873  //bbs[proci][0] = boundBox::invertedBox;
874  bbs[proci][0].min() = point( VGREAT, VGREAT, VGREAT);
875  bbs[proci][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
876  }
877 
878  forAll(s, triI)
879  {
880  point& bbMin = bbs[distribution[triI]][0].min();
881  point& bbMax = bbs[distribution[triI]][0].max();
882 
883  const triSurface::FaceType& f = s[triI];
884  forAll(f, fp)
885  {
886  const point& pt = s.points()[f[fp]];
887  bbMin = ::Foam::min(bbMin, pt);
888  bbMax = ::Foam::max(bbMax, pt);
889  }
890  }
891 
892  // Now combine for all processors and convert to correct format.
893  forAll(bbs, proci)
894  {
895  forAll(bbs[proci], i)
896  {
897  reduce(bbs[proci][i].min(), minOp<point>());
898  reduce(bbs[proci][i].max(), maxOp<point>());
899  }
900  }
901  return bbs;
902 }
903 
904 
905 // Does any part of triangle overlap bb.
906 bool Foam::distributedTriSurfaceMesh::overlaps
907 (
908  const List<treeBoundBox>& bbs,
909  const point& p0,
910  const point& p1,
911  const point& p2
912 )
913 {
914  forAll(bbs, bbI)
915  {
916  const treeBoundBox& bb = bbs[bbI];
917 
918  treeBoundBox triBb(p0, p0);
919  triBb.min() = min(triBb.min(), p1);
920  triBb.min() = min(triBb.min(), p2);
921 
922  triBb.max() = max(triBb.max(), p1);
923  triBb.max() = max(triBb.max(), p2);
924 
925  // Exact test of triangle intersecting bb
926 
927  // Quick rejection. If whole bounding box of tri is outside cubeBb then
928  // there will be no intersection.
929  if (bb.overlaps(triBb))
930  {
931  // Check if one or more triangle point inside
932  if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
933  {
934  // One or more points inside
935  return true;
936  }
937 
938  // Now we have the difficult case: all points are outside but
939  // connecting edges might go through cube. Use fast intersection
940  // of bounding box.
941  bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
942 
943  if (intersect)
944  {
945  return true;
946  }
947  }
948  }
949  return false;
950 }
951 
952 
953 void Foam::distributedTriSurfaceMesh::subsetMeshMap
954 (
955  const triSurface& s,
956  const boolList& include,
957  const label nIncluded,
958  labelList& newToOldPoints,
959  labelList& oldToNewPoints,
960  labelList& newToOldFaces
961 )
962 {
963  newToOldFaces.setSize(nIncluded);
964  newToOldPoints.setSize(s.points().size());
965  oldToNewPoints.setSize(s.points().size());
966  oldToNewPoints = -1;
967  {
968  label facei = 0;
969  label pointi = 0;
970 
971  forAll(include, oldFacei)
972  {
973  if (include[oldFacei])
974  {
975  // Store new faces compact
976  newToOldFaces[facei++] = oldFacei;
977 
978  // Renumber labels for face
979  const triSurface::FaceType& f = s[oldFacei];
980 
981  forAll(f, fp)
982  {
983  label oldPointi = f[fp];
984 
985  if (oldToNewPoints[oldPointi] == -1)
986  {
987  oldToNewPoints[oldPointi] = pointi;
988  newToOldPoints[pointi++] = oldPointi;
989  }
990  }
991  }
992  }
993  newToOldPoints.setSize(pointi);
994  }
995 }
996 
997 
998 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
999 (
1000  const triSurface& s,
1001  const labelList& newToOldPoints,
1002  const labelList& oldToNewPoints,
1003  const labelList& newToOldFaces
1004 )
1005 {
1006  // Extract points
1007  pointField newPoints(newToOldPoints.size());
1008  forAll(newToOldPoints, i)
1009  {
1010  newPoints[i] = s.points()[newToOldPoints[i]];
1011  }
1012  // Extract faces
1013  List<labelledTri> newTriangles(newToOldFaces.size());
1014 
1015  forAll(newToOldFaces, i)
1016  {
1017  // Get old vertex labels
1018  const labelledTri& tri = s[newToOldFaces[i]];
1019 
1020  newTriangles[i][0] = oldToNewPoints[tri[0]];
1021  newTriangles[i][1] = oldToNewPoints[tri[1]];
1022  newTriangles[i][2] = oldToNewPoints[tri[2]];
1023  newTriangles[i].region() = tri.region();
1024  }
1025 
1026  // Reuse storage.
1027  return triSurface(newTriangles, s.patches(), newPoints, true);
1028 }
1029 
1030 
1031 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
1032 (
1033  const triSurface& s,
1034  const boolList& include,
1035  labelList& newToOldPoints,
1036  labelList& newToOldFaces
1037 )
1038 {
1039  label n = 0;
1040 
1041  forAll(include, i)
1042  {
1043  if (include[i])
1044  {
1045  n++;
1046  }
1047  }
1048 
1049  labelList oldToNewPoints;
1050  subsetMeshMap
1051  (
1052  s,
1053  include,
1054  n,
1055  newToOldPoints,
1056  oldToNewPoints,
1057  newToOldFaces
1058  );
1059 
1060  return subsetMesh
1061  (
1062  s,
1063  newToOldPoints,
1064  oldToNewPoints,
1065  newToOldFaces
1066  );
1067 }
1068 
1069 
1070 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
1071 (
1072  const triSurface& s,
1073  const labelList& newToOldFaces,
1074  labelList& newToOldPoints
1075 )
1076 {
1077  const boolList include
1078  (
1079  createWithValues<boolList>
1080  (
1081  s.size(),
1082  false,
1083  newToOldFaces,
1084  true
1085  )
1086  );
1087 
1088  newToOldPoints.setSize(s.points().size());
1089  labelList oldToNewPoints(s.points().size(), -1);
1090  {
1091  label pointi = 0;
1092 
1093  forAll(include, oldFacei)
1094  {
1095  if (include[oldFacei])
1096  {
1097  // Renumber labels for face
1098  const triSurface::FaceType& f = s[oldFacei];
1099 
1100  forAll(f, fp)
1101  {
1102  label oldPointi = f[fp];
1103 
1104  if (oldToNewPoints[oldPointi] == -1)
1105  {
1106  oldToNewPoints[oldPointi] = pointi;
1107  newToOldPoints[pointi++] = oldPointi;
1108  }
1109  }
1110  }
1111  }
1112  newToOldPoints.setSize(pointi);
1113  }
1114 
1115  return subsetMesh
1116  (
1117  s,
1118  newToOldPoints,
1119  oldToNewPoints,
1120  newToOldFaces
1121  );
1122 }
1123 
1124 
1125 Foam::label Foam::distributedTriSurfaceMesh::findTriangle
1126 (
1127  const List<labelledTri>& allFaces,
1128  const labelListList& allPointFaces,
1129  const labelledTri& otherF
1130 )
1131 {
1132  // allFaces connected to otherF[0]
1133  const labelList& pFaces = allPointFaces[otherF[0]];
1134 
1135  forAll(pFaces, i)
1136  {
1137  const labelledTri& f = allFaces[pFaces[i]];
1138 
1139  if (f.region() == otherF.region())
1140  {
1141  // Find index of otherF[0]
1142  label fp0 = findIndex(f, otherF[0]);
1143  // Check rest of triangle in same order
1144  label fp1 = f.fcIndex(fp0);
1145  label fp2 = f.fcIndex(fp1);
1146 
1147  if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
1148  {
1149  return pFaces[i];
1150  }
1151  }
1152  }
1153  return -1;
1154 }
1155 
1156 
1157 // Merge into allSurf.
1158 void Foam::distributedTriSurfaceMesh::merge
1159 (
1160  const scalar mergeDist,
1161  const List<labelledTri>& subTris,
1162  const pointField& subPoints,
1163 
1164  List<labelledTri>& allTris,
1165  pointField& allPoints,
1166 
1167  labelList& faceConstructMap,
1168  labelList& pointConstructMap
1169 )
1170 {
1171  labelList subToAll;
1172  matchPoints
1173  (
1174  subPoints,
1175  allPoints,
1176  scalarField(subPoints.size(), mergeDist), // match distance
1177  false, // verbose
1178  pointConstructMap
1179  );
1180 
1181  label nOldAllPoints = allPoints.size();
1182 
1183 
1184  // Add all unmatched points
1185  // ~~~~~~~~~~~~~~~~~~~~~~~~
1186 
1187  label allPointi = nOldAllPoints;
1188  forAll(pointConstructMap, pointi)
1189  {
1190  if (pointConstructMap[pointi] == -1)
1191  {
1192  pointConstructMap[pointi] = allPointi++;
1193  }
1194  }
1195 
1196  if (allPointi > nOldAllPoints)
1197  {
1198  allPoints.setSize(allPointi);
1199 
1200  forAll(pointConstructMap, pointi)
1201  {
1202  if (pointConstructMap[pointi] >= nOldAllPoints)
1203  {
1204  allPoints[pointConstructMap[pointi]] = subPoints[pointi];
1205  }
1206  }
1207  }
1208 
1209 
1210  // To check whether triangles are same we use pointFaces.
1211  labelListList allPointFaces;
1212  invertManyToMany(nOldAllPoints, allTris, allPointFaces);
1213 
1214 
1215  // Add all unmatched triangles
1216  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1217 
1218  label allTriI = allTris.size();
1219  allTris.setSize(allTriI + subTris.size());
1220 
1221  faceConstructMap.setSize(subTris.size());
1222 
1223  forAll(subTris, triI)
1224  {
1225  const labelledTri& subTri = subTris[triI];
1226 
1227  // Get triangle in new numbering
1228  labelledTri mappedTri
1229  (
1230  pointConstructMap[subTri[0]],
1231  pointConstructMap[subTri[1]],
1232  pointConstructMap[subTri[2]],
1233  subTri.region()
1234  );
1235 
1236 
1237  // Check if all points were already in surface
1238  bool fullMatch = true;
1239 
1240  forAll(mappedTri, fp)
1241  {
1242  if (mappedTri[fp] >= nOldAllPoints)
1243  {
1244  fullMatch = false;
1245  break;
1246  }
1247  }
1248 
1249  if (fullMatch)
1250  {
1251  // All three points are mapped to old points. See if same
1252  // triangle.
1253  label i = findTriangle
1254  (
1255  allTris,
1256  allPointFaces,
1257  mappedTri
1258  );
1259 
1260  if (i == -1)
1261  {
1262  // Add
1263  faceConstructMap[triI] = allTriI;
1264  allTris[allTriI] = mappedTri;
1265  allTriI++;
1266  }
1267  else
1268  {
1269  faceConstructMap[triI] = i;
1270  }
1271  }
1272  else
1273  {
1274  // Add
1275  faceConstructMap[triI] = allTriI;
1276  allTris[allTriI] = mappedTri;
1277  allTriI++;
1278  }
1279  }
1280  allTris.setSize(allTriI);
1281 }
1282 
1283 
1284 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1285 
1286 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
1288  const IOobject& io,
1289  const triSurface& s,
1290  const dictionary& dict
1291 )
1292 :
1293  triSurfaceMesh(io, s),
1294  dict_
1295  (
1296  IOobject
1297  (
1298  searchableSurface::name() + "Dict",
1305  ),
1306  dict
1307  )
1308 {
1309  read();
1310 
1311  reduce(bounds().min(), minOp<point>());
1312  reduce(bounds().max(), maxOp<point>());
1313 
1314  if (debug)
1315  {
1316  InfoInFunction << "Constructed from triSurface:" << endl;
1317  writeStats(Info);
1318 
1319  labelList nTris(Pstream::nProcs());
1320  nTris[Pstream::myProcNo()] = triSurface::size();
1321  Pstream::gatherList(nTris);
1322  Pstream::scatterList(nTris);
1323 
1324  Info<< endl<< "\tproc\ttris\tbb" << endl;
1325  forAll(nTris, proci)
1326  {
1327  Info<< '\t' << proci << '\t' << nTris[proci]
1328  << '\t' << procBb_[proci] << endl;
1329  }
1330  Info<< endl;
1331  }
1332 }
1333 
1334 
1335 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
1336 :
1338  (
1339  IOobject
1340  (
1341  io.name(),
1342  io.time().findInstance(io.local(), word::null),
1343  io.local(),
1344  io.db(),
1345  io.readOpt(),
1346  io.writeOpt(),
1347  io.registerObject()
1348  ),
1349  false
1350  ),
1351  dict_
1352  (
1353  IOobject
1354  (
1355  searchableSurface::name() + "Dict",
1356  searchableSurface::instance(),
1357  searchableSurface::local(),
1358  searchableSurface::db(),
1359  searchableSurface::readOpt(),
1360  searchableSurface::writeOpt(),
1361  searchableSurface::registerObject()
1362  )
1363  )
1364 {
1365  if
1366  (
1367  Pstream::parRun()
1368  && (
1369  dict_.readOpt() == IOobject::MUST_READ
1371  )
1372  && (
1375  )
1376  )
1377  {
1379  << " using 'timeStampMaster' or 'inotifyMaster.'\n"
1380  << " Modify the entry fileModificationChecking\n"
1381  << " in the etc/controlDict.\n"
1382  << " Use 'timeStamp' instead."
1383  << exit(FatalError);
1384  }
1385 
1386  read();
1387 
1388  reduce(bounds().min(), minOp<point>());
1389  reduce(bounds().max(), maxOp<point>());
1390 
1391  if (debug)
1392  {
1393  InfoInFunction << "Read distributedTriSurface from " << io.objectPath()
1394  << ':' << endl;
1395  writeStats(Info);
1396 
1397  labelList nTris(Pstream::nProcs());
1398  nTris[Pstream::myProcNo()] = triSurface::size();
1399  Pstream::gatherList(nTris);
1400  Pstream::scatterList(nTris);
1401 
1402  Info<< endl<< "\tproc\ttris\tbb" << endl;
1403  forAll(nTris, proci)
1404  {
1405  Info<< '\t' << proci << '\t' << nTris[proci]
1406  << '\t' << procBb_[proci] << endl;
1407  }
1408  Info<< endl;
1409  }
1410 }
1411 
1412 
1413 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
1415  const IOobject& io,
1416  const dictionary& dict
1417 )
1418 :
1419  //triSurfaceMesh(io, dict),
1421  (
1422  IOobject
1423  (
1424  io.name(),
1425  io.time().findInstance(io.local(), word::null),
1426  io.local(),
1427  io.db(),
1428  io.readOpt(),
1429  io.writeOpt(),
1430  io.registerObject()
1431  ),
1432  dict,
1433  false
1434  ),
1435  dict_
1436  (
1437  IOobject
1438  (
1439  searchableSurface::name() + "Dict",
1446  )
1447  )
1448 {
1449  if
1450  (
1451  Pstream::parRun()
1452  && (
1453  dict_.readOpt() == IOobject::MUST_READ
1455  )
1456  && (
1459  )
1460  )
1461  {
1463  << " using 'timeStampMaster' or 'inotifyMaster.'\n"
1464  << " Modify the entry fileModificationChecking\n"
1465  << " in the etc/controlDict.\n"
1466  << " Use 'timeStamp' instead."
1467  << exit(FatalError);
1468  }
1469 
1470  read();
1471 
1472  reduce(bounds().min(), minOp<point>());
1473  reduce(bounds().max(), maxOp<point>());
1474 
1475  if (debug)
1476  {
1477  InfoInFunction << "Read distributedTriSurface from " << io.objectPath()
1478  << " and dictionary:" << endl;
1479  writeStats(Info);
1480 
1481  labelList nTris(Pstream::nProcs());
1482  nTris[Pstream::myProcNo()] = triSurface::size();
1483  Pstream::gatherList(nTris);
1484  Pstream::scatterList(nTris);
1485 
1486  Info<< endl<< "\tproc\ttris\tbb" << endl;
1487  forAll(nTris, proci)
1488  {
1489  Info<< '\t' << proci << '\t' << nTris[proci]
1490  << '\t' << procBb_[proci] << endl;
1491  }
1492  Info<< endl;
1493  }
1494 }
1495 
1496 
1497 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1498 
1500 {
1501  clearOut();
1502 }
1503 
1504 
1506 {
1507  globalTris_.clear();
1509 }
1510 
1511 
1512 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1513 
1515 {
1516  if (!globalTris_.valid())
1517  {
1518  globalTris_.reset(new globalIndex(triSurface::size()));
1519  }
1520  return globalTris_;
1521 }
1522 
1523 
1526  const pointField& samples,
1527  const scalarField& nearestDistSqr,
1528  List<pointIndexHit>& info
1529 ) const
1530 {
1531  const indexedOctree<treeDataTriSurface>& octree = tree();
1532 
1533  // Important:force synchronised construction of indexing
1534  const globalIndex& triIndexer = globalTris();
1535 
1536 
1537  // Initialise
1538  // ~~~~~~~~~~
1539 
1540  info.setSize(samples.size());
1541  forAll(info, i)
1542  {
1543  info[i].setMiss();
1544  }
1545 
1546 
1547 
1548  // Do any local queries
1549  // ~~~~~~~~~~~~~~~~~~~~
1550 
1551  label nLocal = 0;
1552 
1553  {
1554  // Work array - whether processor bb overlaps the bounding sphere.
1555  boolList procBbOverlaps(Pstream::nProcs());
1556 
1557  forAll(samples, i)
1558  {
1559  // Find the processor this sample+radius overlaps.
1560  label nProcs = calcOverlappingProcs
1561  (
1562  samples[i],
1563  nearestDistSqr[i],
1564  procBbOverlaps
1565  );
1566 
1567  // Overlaps local processor?
1568  if (procBbOverlaps[Pstream::myProcNo()])
1569  {
1570  info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
1571  if (info[i].hit())
1572  {
1573  info[i].setIndex(triIndexer.toGlobal(info[i].index()));
1574  }
1575  if (nProcs == 1)
1576  {
1577  // Fully local
1578  nLocal++;
1579  }
1580  }
1581  }
1582  }
1583 
1584 
1585  if
1586  (
1587  Pstream::parRun()
1588  && (
1589  returnReduce(nLocal, sumOp<label>())
1590  < returnReduce(samples.size(), sumOp<label>())
1591  )
1592  )
1593  {
1594  // Not all can be resolved locally. Build queries and map, send over
1595  // queries, do intersections, send back and merge.
1596 
1597  // Calculate queries and exchange map
1598  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1599 
1600  pointField allCentres;
1601  scalarField allRadiusSqr;
1602  labelList allSegmentMap;
1603  autoPtr<mapDistribute> mapPtr
1604  (
1605  calcLocalQueries
1606  (
1607  samples,
1608  nearestDistSqr,
1609 
1610  allCentres,
1611  allRadiusSqr,
1612  allSegmentMap
1613  )
1614  );
1615  const mapDistribute& map = mapPtr();
1616 
1617 
1618  // swap samples to local processor
1619  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1620 
1621  map.distribute(allCentres);
1622  map.distribute(allRadiusSqr);
1623 
1624 
1625  // Do my tests
1626  // ~~~~~~~~~~~
1627 
1628  List<pointIndexHit> allInfo(allCentres.size());
1629  forAll(allInfo, i)
1630  {
1631  allInfo[i] = octree.findNearest
1632  (
1633  allCentres[i],
1634  allRadiusSqr[i]
1635  );
1636  if (allInfo[i].hit())
1637  {
1638  allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
1639  }
1640  }
1641 
1642 
1643  // Send back results
1644  // ~~~~~~~~~~~~~~~~~
1645 
1646  map.reverseDistribute(allSegmentMap.size(), allInfo);
1647 
1648 
1649  // Extract information
1650  // ~~~~~~~~~~~~~~~~~~~
1651 
1652  forAll(allInfo, i)
1653  {
1654  if (allInfo[i].hit())
1655  {
1656  label pointi = allSegmentMap[i];
1657 
1658  if (!info[pointi].hit())
1659  {
1660  // No intersection yet so take this one
1661  info[pointi] = allInfo[i];
1662  }
1663  else
1664  {
1665  // Nearest intersection
1666  if
1667  (
1668  magSqr(allInfo[i].hitPoint()-samples[pointi])
1669  < magSqr(info[pointi].hitPoint()-samples[pointi])
1670  )
1671  {
1672  info[pointi] = allInfo[i];
1673  }
1674  }
1675  }
1676  }
1677  }
1678 }
1679 
1680 
1681 void Foam::distributedTriSurfaceMesh::findLine
1683  const pointField& start,
1684  const pointField& end,
1685  List<pointIndexHit>& info
1686 ) const
1687 {
1688  findLine
1689  (
1690  true, // nearestIntersection
1691  start,
1692  end,
1693  info
1694  );
1695 }
1696 
1697 
1700  const pointField& start,
1701  const pointField& end,
1702  List<pointIndexHit>& info
1703 ) const
1704 {
1705  findLine
1706  (
1707  true, // nearestIntersection
1708  start,
1709  end,
1710  info
1711  );
1712 }
1713 
1714 
1717  const pointField& start,
1718  const pointField& end,
1719  List<List<pointIndexHit>>& info
1720 ) const
1721 {
1722  // Reuse fineLine. We could modify all of findLine to do multiple
1723  // intersections but this would mean a lot of data transferred so
1724  // for now we just find nearest intersection and retest from that
1725  // intersection onwards.
1726 
1727  // Work array.
1728  List<pointIndexHit> hitInfo(start.size());
1729 
1730  findLine
1731  (
1732  true, // nearestIntersection
1733  start,
1734  end,
1735  hitInfo
1736  );
1737 
1738  // Tolerances:
1739  // To find all intersections we add a small vector to the last intersection
1740  // This is chosen such that
1741  // - it is significant (SMALL is smallest representative relative tolerance;
1742  // we need something bigger since we're doing calculations)
1743  // - if the start-end vector is zero we still progress
1744  const vectorField dirVec(end-start);
1745  const scalarField magSqrDirVec(magSqr(dirVec));
1746  const vectorField smallVec
1747  (
1748  ROOTSMALL*dirVec
1749  + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
1750  );
1751 
1752  // Copy to input and compact any hits
1753  labelList pointMap(start.size());
1754  pointField e0(start.size());
1755  pointField e1(start.size());
1756  label compactI = 0;
1757 
1758  info.setSize(hitInfo.size());
1759  forAll(hitInfo, pointi)
1760  {
1761  if (hitInfo[pointi].hit())
1762  {
1763  info[pointi].setSize(1);
1764  info[pointi][0] = hitInfo[pointi];
1765 
1766  point pt = hitInfo[pointi].hitPoint() + smallVec[pointi];
1767 
1768  if (((pt-start[pointi])&dirVec[pointi]) <= magSqrDirVec[pointi])
1769  {
1770  e0[compactI] = pt;
1771  e1[compactI] = end[pointi];
1772  pointMap[compactI] = pointi;
1773  compactI++;
1774  }
1775  }
1776  else
1777  {
1778  info[pointi].clear();
1779  }
1780  }
1781 
1782  e0.setSize(compactI);
1783  e1.setSize(compactI);
1784  pointMap.setSize(compactI);
1785 
1786  while (returnReduce(e0.size(), sumOp<label>()) > 0)
1787  {
1788  findLine
1789  (
1790  true, // nearestIntersection
1791  e0,
1792  e1,
1793  hitInfo
1794  );
1795 
1796 
1797  // Extract
1798  label compactI = 0;
1799  forAll(hitInfo, i)
1800  {
1801  if (hitInfo[i].hit())
1802  {
1803  label pointi = pointMap[i];
1804 
1805  label sz = info[pointi].size();
1806  info[pointi].setSize(sz+1);
1807  info[pointi][sz] = hitInfo[i];
1808 
1809  point pt = hitInfo[i].hitPoint() + smallVec[pointi];
1810 
1811  if (((pt-start[pointi])&dirVec[pointi]) <= magSqrDirVec[pointi])
1812  {
1813  e0[compactI] = pt;
1814  e1[compactI] = end[pointi];
1815  pointMap[compactI] = pointi;
1816  compactI++;
1817  }
1818  }
1819  }
1820 
1821  // Trim
1822  e0.setSize(compactI);
1823  e1.setSize(compactI);
1824  pointMap.setSize(compactI);
1825  }
1826 }
1827 
1828 
1831  const List<pointIndexHit>& info,
1832  labelList& region
1833 ) const
1834 {
1835  if (!Pstream::parRun())
1836  {
1837  region.setSize(info.size());
1838  forAll(info, i)
1839  {
1840  if (info[i].hit())
1841  {
1842  region[i] = triSurface::operator[](info[i].index()).region();
1843  }
1844  else
1845  {
1846  region[i] = -1;
1847  }
1848  }
1849  return;
1850  }
1851 
1852 
1853  // Get query data (= local index of triangle)
1854  // ~~~~~~~~~~~~~~
1855 
1856  labelList triangleIndex(info.size());
1857  autoPtr<mapDistribute> mapPtr
1858  (
1859  calcLocalQueries
1860  (
1861  info,
1862  triangleIndex
1863  )
1864  );
1865  const mapDistribute& map = mapPtr();
1866 
1867 
1868  // Do my tests
1869  // ~~~~~~~~~~~
1870 
1871  const triSurface& s = static_cast<const triSurface&>(*this);
1872 
1873  region.setSize(triangleIndex.size());
1874 
1875  forAll(triangleIndex, i)
1876  {
1877  label triI = triangleIndex[i];
1878  region[i] = s[triI].region();
1879  }
1880 
1881 
1882  // Send back results
1883  // ~~~~~~~~~~~~~~~~~
1884 
1885  map.reverseDistribute(info.size(), region);
1886 }
1887 
1888 
1891  const List<pointIndexHit>& info,
1892  vectorField& normal
1893 ) const
1894 {
1895  if (!Pstream::parRun())
1896  {
1897  triSurfaceMesh::getNormal(info, normal);
1898  return;
1899  }
1900 
1901 
1902  // Get query data (= local index of triangle)
1903  // ~~~~~~~~~~~~~~
1904 
1905  labelList triangleIndex(info.size());
1906  autoPtr<mapDistribute> mapPtr
1907  (
1908  calcLocalQueries
1909  (
1910  info,
1911  triangleIndex
1912  )
1913  );
1914  const mapDistribute& map = mapPtr();
1915 
1916 
1917  // Do my tests
1918  // ~~~~~~~~~~~
1919 
1920  const triSurface& s = static_cast<const triSurface&>(*this);
1921 
1922  normal.setSize(triangleIndex.size());
1923 
1924  forAll(triangleIndex, i)
1925  {
1926  label triI = triangleIndex[i];
1927  normal[i] = s[triI].normal(s.points());
1928  normal[i] /= mag(normal[i]) + VSMALL;
1929  }
1930 
1931 
1932  // Send back results
1933  // ~~~~~~~~~~~~~~~~~
1934 
1935  map.reverseDistribute(info.size(), normal);
1936 }
1937 
1938 
1941  const List<pointIndexHit>& info,
1942  labelList& values
1943 ) const
1944 {
1945  if (!Pstream::parRun())
1946  {
1947  triSurfaceMesh::getField(info, values);
1948  return;
1949  }
1950 
1951  if (foundObject<triSurfaceLabelField>("values"))
1952  {
1953  const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
1954  (
1955  "values"
1956  );
1957 
1958 
1959  // Get query data (= local index of triangle)
1960  // ~~~~~~~~~~~~~~
1961 
1962  labelList triangleIndex(info.size());
1963  autoPtr<mapDistribute> mapPtr
1964  (
1965  calcLocalQueries
1966  (
1967  info,
1968  triangleIndex
1969  )
1970  );
1971  const mapDistribute& map = mapPtr();
1972 
1973 
1974  // Do my tests
1975  // ~~~~~~~~~~~
1976 
1977  values.setSize(triangleIndex.size());
1978 
1979  forAll(triangleIndex, i)
1980  {
1981  label triI = triangleIndex[i];
1982  values[i] = fld[triI];
1983  }
1984 
1985 
1986  // Send back results
1987  // ~~~~~~~~~~~~~~~~~
1988 
1989  map.reverseDistribute(info.size(), values);
1990  }
1991 }
1992 
1993 
1996  const pointField& points,
1997  List<volumeType>& volType
1998 ) const
1999 {
2001  << "Volume type not supported for distributed surfaces."
2002  << exit(FatalError);
2003 }
2004 
2005 
2006 // Subset the part of surface that is overlapping bb.
2009  const triSurface& s,
2010  const List<treeBoundBox>& bbs,
2011 
2012  labelList& subPointMap,
2013  labelList& subFaceMap
2014 )
2015 {
2016  // Determine what triangles to keep.
2017  boolList includedFace(s.size(), false);
2018 
2019  // Create a slightly larger bounding box.
2020  List<treeBoundBox> bbsX(bbs.size());
2021  const scalar eps = 1.0e-4;
2022  forAll(bbs, i)
2023  {
2024  const point mid = 0.5*(bbs[i].min() + bbs[i].max());
2025  const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
2026 
2027  bbsX[i].min() = mid - halfSpan;
2028  bbsX[i].max() = mid + halfSpan;
2029  }
2030 
2031  forAll(s, triI)
2032  {
2033  const labelledTri& f = s[triI];
2034  const point& p0 = s.points()[f[0]];
2035  const point& p1 = s.points()[f[1]];
2036  const point& p2 = s.points()[f[2]];
2037 
2038  if (overlaps(bbsX, p0, p1, p2))
2039  {
2040  includedFace[triI] = true;
2041  }
2042  }
2043 
2044  return subsetMesh(s, includedFace, subPointMap, subFaceMap);
2045 }
2046 
2047 
2050  const List<treeBoundBox>& bbs,
2051  const bool keepNonLocal,
2052  autoPtr<mapDistribute>& faceMap,
2053  autoPtr<mapDistribute>& pointMap
2054 )
2055 {
2056  // Get bbs of all domains
2057  // ~~~~~~~~~~~~~~~~~~~~~~
2058 
2059  {
2061 
2062  switch(distType_)
2063  {
2064  case FOLLOW:
2065  newProcBb[Pstream::myProcNo()].setSize(bbs.size());
2066  forAll(bbs, i)
2067  {
2068  newProcBb[Pstream::myProcNo()][i] = bbs[i];
2069  }
2070  Pstream::gatherList(newProcBb);
2071  Pstream::scatterList(newProcBb);
2072  break;
2073 
2074  case INDEPENDENT:
2075  newProcBb = independentlyDistributedBbs(*this);
2076  break;
2077 
2078  case FROZEN:
2079  return;
2080  break;
2081 
2082  default:
2084  << "Unsupported distribution type." << exit(FatalError);
2085  break;
2086  }
2087 
2088  if (newProcBb == procBb_)
2089  {
2090  return;
2091  }
2092  else
2093  {
2094  procBb_.transfer(newProcBb);
2095  dict_.set("bounds", procBb_[Pstream::myProcNo()]);
2096  }
2097  }
2098 
2099 
2100  // Debug information
2101  if (debug)
2102  {
2103  labelList nTris(Pstream::nProcs());
2104  nTris[Pstream::myProcNo()] = triSurface::size();
2105  Pstream::gatherList(nTris);
2106  Pstream::scatterList(nTris);
2107 
2109  << "before distribution:" << endl << "\tproc\ttris" << endl;
2110 
2111  forAll(nTris, proci)
2112  {
2113  Info<< '\t' << proci << '\t' << nTris[proci] << endl;
2114  }
2115  Info<< endl;
2116  }
2117 
2118 
2119  // Use procBbs to determine which faces go where
2120  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2121 
2122  labelListList faceSendMap(Pstream::nProcs());
2123  labelListList pointSendMap(Pstream::nProcs());
2124 
2125  forAll(procBb_, proci)
2126  {
2128  (
2129  *this,
2130  procBb_[proci],
2131  pointSendMap[proci],
2132  faceSendMap[proci]
2133  );
2134 
2135  if (debug)
2136  {
2137  //Pout<< "Overlapping with proc " << proci
2138  // << " faces:" << faceSendMap[proci].size()
2139  // << " points:" << pointSendMap[proci].size() << endl << endl;
2140  }
2141  }
2142 
2143  if (keepNonLocal)
2144  {
2145  // Include in faceSendMap/pointSendMap the triangles that are
2146  // not mapped to any processor so they stay local.
2147 
2148  const triSurface& s = static_cast<const triSurface&>(*this);
2149 
2150  boolList includedFace(s.size(), true);
2151 
2152  forAll(faceSendMap, proci)
2153  {
2154  if (proci != Pstream::myProcNo())
2155  {
2156  forAll(faceSendMap[proci], i)
2157  {
2158  includedFace[faceSendMap[proci][i]] = false;
2159  }
2160  }
2161  }
2162 
2163  // Combine includedFace (all triangles that are not on any neighbour)
2164  // with overlap.
2165 
2166  forAll(faceSendMap[Pstream::myProcNo()], i)
2167  {
2168  includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
2169  }
2170 
2171  subsetMesh
2172  (
2173  s,
2174  includedFace,
2175  pointSendMap[Pstream::myProcNo()],
2176  faceSendMap[Pstream::myProcNo()]
2177  );
2178  }
2179 
2180 
2181  // Send over how many faces/points I need to receive
2182  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2183 
2184  labelListList faceSendSizes(Pstream::nProcs());
2185  faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
2186  forAll(faceSendMap, proci)
2187  {
2188  faceSendSizes[Pstream::myProcNo()][proci] = faceSendMap[proci].size();
2189  }
2190  Pstream::gatherList(faceSendSizes);
2191  Pstream::scatterList(faceSendSizes);
2192 
2193 
2194 
2195  // Exchange surfaces
2196  // ~~~~~~~~~~~~~~~~~
2197 
2198  // Storage for resulting surface
2199  List<labelledTri> allTris;
2201 
2202  labelListList faceConstructMap(Pstream::nProcs());
2203  labelListList pointConstructMap(Pstream::nProcs());
2204 
2205 
2206  // My own surface first
2207  // ~~~~~~~~~~~~~~~~~~~~
2208 
2209  {
2210  labelList pointMap;
2211  triSurface subSurface
2212  (
2213  subsetMesh
2214  (
2215  *this,
2216  faceSendMap[Pstream::myProcNo()],
2217  pointMap
2218  )
2219  );
2220 
2221  allTris = subSurface;
2222  allPoints = subSurface.points();
2223 
2224  faceConstructMap[Pstream::myProcNo()] = identity
2225  (
2226  faceSendMap[Pstream::myProcNo()].size()
2227  );
2228  pointConstructMap[Pstream::myProcNo()] = identity
2229  (
2230  pointSendMap[Pstream::myProcNo()].size()
2231  );
2232  }
2233 
2234 
2235 
2236  // Send all
2237  // ~~~~~~~~
2238 
2239  forAll(faceSendSizes, proci)
2240  {
2241  if (proci != Pstream::myProcNo())
2242  {
2243  if (faceSendSizes[Pstream::myProcNo()][proci] > 0)
2244  {
2246 
2247  labelList pointMap;
2248  triSurface subSurface
2249  (
2250  subsetMesh
2251  (
2252  *this,
2253  faceSendMap[proci],
2254  pointMap
2255  )
2256  );
2257 
2258  //if (debug)
2259  //{
2260  // Pout<< "Sending to " << proci
2261  // << " faces:" << faceSendMap[proci].size()
2262  // << " points:" << subSurface.points().size() << endl
2263  // << endl;
2264  //}
2265 
2266  str << subSurface;
2267  }
2268  }
2269  }
2270 
2271 
2272  // Receive and merge all
2273  // ~~~~~~~~~~~~~~~~~~~~~
2274 
2275  forAll(faceSendSizes, proci)
2276  {
2277  if (proci != Pstream::myProcNo())
2278  {
2279  if (faceSendSizes[proci][Pstream::myProcNo()] > 0)
2280  {
2282 
2283  // Receive
2284  triSurface subSurface(str);
2285 
2286  //if (debug)
2287  //{
2288  // Pout<< "Received from " << proci
2289  // << " faces:" << subSurface.size()
2290  // << " points:" << subSurface.points().size() << endl
2291  // << endl;
2292  //}
2293 
2294  // Merge into allSurf
2295  merge
2296  (
2297  mergeDist_,
2298  subSurface,
2299  subSurface.points(),
2300 
2301  allTris,
2302  allPoints,
2303  faceConstructMap[proci],
2304  pointConstructMap[proci]
2305  );
2306 
2307  //if (debug)
2308  //{
2309  // Pout<< "Current merged surface : faces:" << allTris.size()
2310  // << " points:" << allPoints.size() << endl << endl;
2311  //}
2312  }
2313  }
2314  }
2315 
2316 
2317  faceMap.reset
2318  (
2319  new mapDistribute
2320  (
2321  allTris.size(),
2322  faceSendMap.xfer(),
2323  faceConstructMap.xfer()
2324  )
2325  );
2326  pointMap.reset
2327  (
2328  new mapDistribute
2329  (
2330  allPoints.size(),
2331  pointSendMap.xfer(),
2332  pointConstructMap.xfer()
2333  )
2334  );
2335 
2336  // Construct triSurface. Reuse storage.
2337  triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
2338 
2339  clearOut();
2340 
2341  // Set the bounds() value to the boundBox of the undecomposed surface
2343 
2344  reduce(bounds().min(), minOp<point>());
2345  reduce(bounds().max(), maxOp<point>());
2346 
2347  // Regions stays same
2348  // Volume type stays same.
2349 
2350  distributeFields<label>(faceMap());
2351  distributeFields<scalar>(faceMap());
2352  distributeFields<vector>(faceMap());
2353  distributeFields<sphericalTensor>(faceMap());
2354  distributeFields<symmTensor>(faceMap());
2355  distributeFields<tensor>(faceMap());
2356 
2357  if (debug)
2358  {
2359  labelList nTris(Pstream::nProcs());
2360  nTris[Pstream::myProcNo()] = triSurface::size();
2361  Pstream::gatherList(nTris);
2362  Pstream::scatterList(nTris);
2363 
2365  << "after distribution:" << endl << "\tproc\ttris" << endl;
2366 
2367  forAll(nTris, proci)
2368  {
2369  Info<< '\t' << proci << '\t' << nTris[proci] << endl;
2370  }
2371  Info<< endl;
2372  }
2373 }
2374 
2375 
2381  const bool valid
2382 ) const
2383 {
2384  // Make sure dictionary goes to same directory as surface
2385  const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
2386 
2387  // Copy of triSurfaceMesh::writeObject except for the sorting of
2388  // triangles by region. This is done so we preserve region names,
2389  // even if locally we have zero triangles.
2390  {
2392 
2393  if (!mkDir(fullPath.path()))
2394  {
2395  return false;
2396  }
2397 
2398  // Important: preserve any zero-sized patches
2399  triSurface::write(fullPath, true);
2400 
2401  if (!isFile(fullPath))
2402  {
2403  return false;
2404  }
2405  }
2406 
2407  // Dictionary needs to be written in ascii - binary output not supported.
2408  return dict_.writeObject(IOstream::ASCII, ver, cmp, true);
2409 }
2410 
2411 
2413 {
2414  boundBox bb;
2415  label nPoints;
2416  PatchTools::calcBounds(static_cast<const triSurface&>(*this), bb, nPoints);
2417  reduce(bb.min(), minOp<point>());
2418  reduce(bb.max(), maxOp<point>());
2419 
2420  os << "Triangles : " << returnReduce(triSurface::size(), sumOp<label>())
2421  << endl
2422  << "Vertices : " << returnReduce(nPoints, sumOp<label>()) << endl
2423  << "Bounding Box : " << bb << endl;
2424 }
2425 
2426 
2427 // ************************************************************************* //
virtual void getField(const List< pointIndexHit > &, labelList &) const
WIP. From a set of hits (points and.
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
label nPoints() const
Return number of points supporting patch faces.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
const word & name() const
Return name.
Definition: IOobject.H:291
A class for handling file names.
Definition: fileName.H:69
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
static const Form max
Definition: VectorSpace.H:112
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void clearOut()
Clear storage.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void operator=(const triSurface &)
Definition: triSurface.C:1102
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
T & operator[](const label)
Return element of UList.
Definition: UListI.H:167
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
Xfer< List< T > > xfer()
Transfer contents to the Xfer container.
Definition: ListI.H:177
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
virtual tmp< pointField > points() const
Get the points that define the surface.
const globalIndex & globalTris() const
Triangle indexing (demand driven)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
void setIndex(const label index)
static bool intersectBb(const point &p0, const point &p1, const point &p2, const treeBoundBox &cubeBb)
Does triangle intersect bounding box.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:418
virtual void distribute(const List< treeBoundBox > &, const bool keepNonLocal, autoPtr< mapDistribute > &faceMap, autoPtr< mapDistribute > &pointMap)
Set bounds of surface. Bounds currently set as list of.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
word findInstance(const fileName &dir, const word &name=word::null, const IOobject::readOption rOpt=IOobject::MUST_READ, const word &stopInstance=word::null) const
Return the location of "dir" containing the file "name".
Definition: findInstance.C:82
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Pair< point > segment
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static const NamedEnum< distributionType, 3 > distributionTypeNames_
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:52
Determine correspondence between points. See below.
Macros for easy insertion into run-time selection tables.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.
static const Form min
Definition: VectorSpace.H:113
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
IOoject and searching on triSurface.
Input inter-processor communications stream.
Definition: IPstream.H:50
virtual bool writeObject(IOstream::streamFormat, IOstream::versionNumber, IOstream::compressionType, const bool valid) const
Write using given format, version and compression.
const Point & hitPoint() const
Return hit point.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
label region() const
Return region label.
Definition: labelledTriI.H:68
bool read(const char *, int32_t &)
Definition: int32IO.C:85
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:33
bool intersects(const point &overallStart, const vector &overallVec, const point &start, const point &end, point &pt, direction &ptBits) const
Intersects segment; set point to intersection position and face,.
Definition: treeBoundBox.C:235
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< ' ';}gmvFile<< nl;forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Geometrical domain decomposition.
Definition: geomDecomp.H:47
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
const fileName & local() const
Definition: IOobject.H:396
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
streamFormat
Enumeration for the format of data in the stream.
Definition: IOstream.H:86
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
bool hit() const
Is there a hit.
static void calcBounds(const PrimitivePatch< Face, FaceList, PointField, PointType > &p, boundBox &bb, label &nPoints)
const Field< PointType > & points() const
Return reference to global points.
static const word null
An empty word.
Definition: word.H:77
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
label toLocal(const label i) const
From global to local on current processor.
Definition: globalIndexI.H:117
Triangle with additional region number.
Definition: labelledTri.H:57
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
label whichProcID(const label i) const
Which processor does global come from? Binary search.
Definition: globalIndexI.H:123
Accumulating histogram of values. Specified bin resolution automatic generation of bins...
Definition: distribution.H:58
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: FixedListI.H:135
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit >> &) const
Get all intersections in order from start to end.
static autoPtr< decompositionMethod > New(const dictionary &decompositionDict)
Return a reference to the selected decomposition method.
compressionType
Enumeration for the format of data in the stream.
Definition: IOstream.H:193
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
Definition: treeBoundBox.C:393
bool isFile(const fileName &, const bool checkGzip=true, const bool followLink=true)
Does the name exist as a FILE in the file system?
Definition: POSIX.C:551
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:240
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
defineTypeNameAndDebug(combustionModel, 0)
static fileCheckTypes fileModificationChecking
Type of file modification checking.
Definition: IOobject.H:208
labelList f(nPoints)
Output inter-processor communications stream.
Definition: OPstream.H:50
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
virtual labelList decompose(const pointField &points, const scalarField &pointWeights)=0
Return for every coordinate the wanted processor number.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static triSurface overlappingSurface(const triSurface &, const List< treeBoundBox > &, labelList &subPointMap, labelList &subFaceMap)
Subset the part of surface that is overlapping bounds.
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:297
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
virtual bool writeObject(IOstream::streamFormat fmt, IOstream::versionNumber ver, IOstream::compressionType cmp, const bool valid) const
Write using given format, version and compression.
void reverseDistribute(const label constructSize, List< T > &, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Reverse distribute data using default commsType.
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
const geometricSurfacePatchList & patches() const
Definition: triSurface.H:314
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
writeOption writeOpt() const
Definition: IOobject.H:363
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:394
const fileName & instance() const
Definition: IOobject.H:386
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:400
triSurface()
Construct null.
Definition: triSurface.C:596
Class containing processor-to-processor mapping information.
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
virtual label size() const
Range of local indices that can be returned.
const Time & time() const
Return time.
Definition: IOobject.C:337
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
void invertManyToMany(const label len, const UList< InList > &, List< OutList > &)
Invert many-to-many.
Version number type.
Definition: IOstream.H:96
void set(entry *)
Assign a new entry, overwrite any existing entry.
Definition: dictionary.C:941
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
fileName path() const
Return directory path name (part before last /)
Definition: fileName.C:250
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const boundBox & bounds() const
Return const reference to boundBox.
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:331
Triangulated surface description with patch information.
Definition: triSurface.H:65
bool & registerObject()
Register object created from this IOobject with registry if true.
Definition: IOobject.H:327
readOption readOpt() const
Definition: IOobject.H:353
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
label size() const
Return the number of elements in the UList.
Definition: ListI.H:170
fileName objectPath() const
Return complete path + object name.
Definition: IOobject.H:412
const indexedOctree< treeDataTriSurface > & tree() const
Demand driven construction of the octree.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
Namespace for OpenFOAM.
virtual void getField(const List< pointIndexHit > &, labelList &) const
WIP. From a set of hits (points and.
int system(const std::string &command)
Execute the specified command.
Definition: POSIX.C:1193
void writeStats(Ostream &os) const
Print some stats. Parallel aware version of.
#define InfoInFunction
Report an information message using Foam::Info.