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