00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "distributedTriSurfaceMesh.H"
00027 #include <OpenFOAM/mapDistribute.H>
00028 #include <OpenFOAM/Random.H>
00029 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00030 #include <meshTools/triangleFuncs.H>
00031 #include <OpenFOAM/matchPoints.H>
00032 #include <OpenFOAM/globalIndex.H>
00033 #include <OpenFOAM/Time.H>
00034
00035 #include <OpenFOAM/IFstream.H>
00036 #include <decompositionMethods/decompositionMethod.H>
00037 #include <OpenFOAM/vectorList.H>
00038 #include <OpenFOAM/PackedBoolList.H>
00039
00040
00041
00042 namespace Foam
00043 {
00044 defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
00045 addToRunTimeSelectionTable
00046 (
00047 searchableSurface,
00048 distributedTriSurfaceMesh,
00049 dict
00050 );
00051 }
00052
00053
00054 template<>
00055 const char*
00056 Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>::names[] =
00057 {
00058 "follow",
00059 "independent",
00060 "frozen"
00061 };
00062
00063 const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
00064 Foam::distributedTriSurfaceMesh::distributionTypeNames_;
00065
00066
00067
00068
00069
00070 bool Foam::distributedTriSurfaceMesh::read()
00071 {
00072
00073 procBb_.setSize(Pstream::nProcs());
00074
00075 procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
00076 Pstream::gatherList(procBb_);
00077 Pstream::scatterList(procBb_);
00078
00079
00080 distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
00081
00082
00083 mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
00084
00085 return true;
00086 }
00087
00088
00089
00090 bool Foam::distributedTriSurfaceMesh::isLocal
00091 (
00092 const List<treeBoundBox>& myBbs,
00093 const point& start,
00094 const point& end
00095 )
00096 {
00097 forAll(myBbs, bbI)
00098 {
00099 if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
00100 {
00101 return true;
00102 }
00103 }
00104 return false;
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 void Foam::distributedTriSurfaceMesh::distributeSegment
00171 (
00172 const label segmentI,
00173 const point& start,
00174 const point& end,
00175
00176 DynamicList<segment>& allSegments,
00177 DynamicList<label>& allSegmentMap,
00178 List<DynamicList<label> >& sendMap
00179 ) const
00180 {
00181
00182 point clipPt;
00183
00184
00185
00186 if (isLocal(procBb_[Pstream::myProcNo()], start, end))
00187 {
00188 return;
00189 }
00190
00191
00192
00193
00194
00195 forAll(procBb_, procI)
00196 {
00197 if (procI != Pstream::myProcNo())
00198 {
00199 const List<treeBoundBox>& bbs = procBb_[procI];
00200
00201 if (isLocal(bbs, start, end))
00202 {
00203 sendMap[procI].append(allSegments.size());
00204 allSegmentMap.append(segmentI);
00205 allSegments.append(segment(start, end));
00206 return;
00207 }
00208 }
00209 }
00210
00211
00212
00213
00214 forAll(procBb_, procI)
00215 {
00216 const List<treeBoundBox>& bbs = procBb_[procI];
00217
00218 forAll(bbs, bbI)
00219 {
00220 const treeBoundBox& bb = bbs[bbI];
00221
00222
00223
00224
00225 if (bb.intersects(start, end, clipPt))
00226 {
00227 sendMap[procI].append(allSegments.size());
00228 allSegmentMap.append(segmentI);
00229 allSegments.append(segment(start, end));
00230 }
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246 }
00247 }
00248 }
00249
00250
00251 Foam::autoPtr<Foam::mapDistribute>
00252 Foam::distributedTriSurfaceMesh::distributeSegments
00253 (
00254 const pointField& start,
00255 const pointField& end,
00256
00257 List<segment>& allSegments,
00258 labelList& allSegmentMap
00259 ) const
00260 {
00261
00262
00263
00264 labelListList sendMap(Pstream::nProcs());
00265
00266 {
00267
00268
00269
00270
00271
00272 DynamicList<segment> dynAllSegments(start.size());
00273
00274 DynamicList<label> dynAllSegmentMap(start.size());
00275
00276 List<DynamicList<label> > dynSendMap(Pstream::nProcs());
00277
00278 forAll(start, segmentI)
00279 {
00280 distributeSegment
00281 (
00282 segmentI,
00283 start[segmentI],
00284 end[segmentI],
00285
00286 dynAllSegments,
00287 dynAllSegmentMap,
00288 dynSendMap
00289 );
00290 }
00291
00292
00293 sendMap.setSize(Pstream::nProcs());
00294 forAll(sendMap, procI)
00295 {
00296 dynSendMap[procI].shrink();
00297 sendMap[procI].transfer(dynSendMap[procI]);
00298 }
00299
00300 allSegments.transfer(dynAllSegments.shrink());
00301 allSegmentMap.transfer(dynAllSegmentMap.shrink());
00302 }
00303
00304
00305
00306 labelListList sendSizes(Pstream::nProcs());
00307 sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00308 forAll(sendMap, procI)
00309 {
00310 sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00311 }
00312 Pstream::gatherList(sendSizes);
00313 Pstream::scatterList(sendSizes);
00314
00315
00316
00317 labelListList constructMap(Pstream::nProcs());
00318
00319
00320 constructMap[Pstream::myProcNo()] = identity
00321 (
00322 sendMap[Pstream::myProcNo()].size()
00323 );
00324
00325 label segmentI = constructMap[Pstream::myProcNo()].size();
00326 forAll(constructMap, procI)
00327 {
00328 if (procI != Pstream::myProcNo())
00329 {
00330
00331 label nRecv = sendSizes[procI][Pstream::myProcNo()];
00332 constructMap[procI].setSize(nRecv);
00333
00334 for (label i = 0; i < nRecv; i++)
00335 {
00336 constructMap[procI][i] = segmentI++;
00337 }
00338 }
00339 }
00340
00341 return autoPtr<mapDistribute>
00342 (
00343 new mapDistribute
00344 (
00345 segmentI,
00346 sendMap,
00347 constructMap,
00348 true
00349 )
00350 );
00351 }
00352
00353
00354 void Foam::distributedTriSurfaceMesh::findLine
00355 (
00356 const bool nearestIntersection,
00357 const pointField& start,
00358 const pointField& end,
00359 List<pointIndexHit>& info
00360 ) const
00361 {
00362 const indexedOctree<treeDataTriSurface>& octree = tree();
00363
00364
00365 const globalIndex& triIndexer = globalTris();
00366
00367
00368 info.setSize(start.size());
00369 forAll(info, i)
00370 {
00371 info[i].setMiss();
00372 }
00373
00374
00375
00376
00377
00378 label nLocal = 0;
00379
00380 forAll(start, i)
00381 {
00382 if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
00383 {
00384 if (nearestIntersection)
00385 {
00386 info[i] = octree.findLine(start[i], end[i]);
00387 }
00388 else
00389 {
00390 info[i] = octree.findLineAny(start[i], end[i]);
00391 }
00392
00393 if (info[i].hit())
00394 {
00395 info[i].setIndex(triIndexer.toGlobal(info[i].index()));
00396 }
00397 nLocal++;
00398 }
00399 }
00400
00401
00402 if
00403 (
00404 Pstream::parRun()
00405 && (
00406 returnReduce(nLocal, sumOp<label>())
00407 < returnReduce(start.size(), sumOp<label>())
00408 )
00409 )
00410 {
00411
00412
00413
00414
00415
00416
00417
00418
00419 List<segment> allSegments(start.size());
00420
00421 labelList allSegmentMap(start.size());
00422
00423 const autoPtr<mapDistribute> mapPtr
00424 (
00425 distributeSegments
00426 (
00427 start,
00428 end,
00429 allSegments,
00430 allSegmentMap
00431 )
00432 );
00433 const mapDistribute& map = mapPtr();
00434
00435 label nOldAllSegments = allSegments.size();
00436
00437
00438
00439
00440
00441 map.distribute
00442 (
00443 Pstream::nonBlocking,
00444 List<labelPair>(0),
00445 map.constructSize(),
00446 map.subMap(),
00447 map.constructMap(),
00448 allSegments
00449 );
00450
00451
00452
00453
00454
00455
00456 List<pointIndexHit> intersections(allSegments.size());
00457
00458 forAll(allSegments, i)
00459 {
00460 if (nearestIntersection)
00461 {
00462 intersections[i] = octree.findLine
00463 (
00464 allSegments[i].first(),
00465 allSegments[i].second()
00466 );
00467 }
00468 else
00469 {
00470 intersections[i] = octree.findLineAny
00471 (
00472 allSegments[i].first(),
00473 allSegments[i].second()
00474 );
00475 }
00476
00477
00478 if (intersections[i].hit())
00479 {
00480 intersections[i].setIndex
00481 (
00482 triIndexer.toGlobal(intersections[i].index())
00483 );
00484 }
00485 }
00486
00487
00488
00489
00490
00491 map.distribute
00492 (
00493
00494
00495
00496
00497
00498
00499 Pstream::nonBlocking,
00500 List<labelPair>(0),
00501 nOldAllSegments,
00502 map.constructMap(),
00503 map.subMap(),
00504 intersections
00505 );
00506
00507
00508
00509
00510
00511 forAll(intersections, i)
00512 {
00513 const pointIndexHit& allInfo = intersections[i];
00514 label segmentI = allSegmentMap[i];
00515 pointIndexHit& hitInfo = info[segmentI];
00516
00517 if (allInfo.hit())
00518 {
00519 if (!hitInfo.hit())
00520 {
00521
00522 hitInfo = allInfo;
00523 }
00524 else if (nearestIntersection)
00525 {
00526
00527 if
00528 (
00529 magSqr(allInfo.hitPoint()-start[segmentI])
00530 < magSqr(hitInfo.hitPoint()-start[segmentI])
00531 )
00532 {
00533 hitInfo = allInfo;
00534 }
00535 }
00536 }
00537 }
00538 }
00539 }
00540
00541
00542
00543
00544
00545 Foam::autoPtr<Foam::mapDistribute>
00546 Foam::distributedTriSurfaceMesh::calcLocalQueries
00547 (
00548 const List<pointIndexHit>& info,
00549 labelList& triangleIndex
00550 ) const
00551 {
00552 triangleIndex.setSize(info.size());
00553
00554 const globalIndex& triIndexer = globalTris();
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 labelList nSend(Pstream::nProcs(), 0);
00565
00566 forAll(info, i)
00567 {
00568 if (info[i].hit())
00569 {
00570 label procI = triIndexer.whichProcID(info[i].index());
00571 nSend[procI]++;
00572 }
00573 }
00574
00575
00576 labelListList sendMap(Pstream::nProcs());
00577 forAll(nSend, procI)
00578 {
00579 sendMap[procI].setSize(nSend[procI]);
00580 nSend[procI] = 0;
00581 }
00582
00583
00584 forAll(info, i)
00585 {
00586 if (info[i].hit())
00587 {
00588 label procI = triIndexer.whichProcID(info[i].index());
00589 triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
00590 sendMap[procI][nSend[procI]++] = i;
00591 }
00592 else
00593 {
00594 triangleIndex[i] = -1;
00595 }
00596 }
00597
00598
00599
00600
00601
00602 labelListList sendSizes(Pstream::nProcs());
00603 sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00604 forAll(sendMap, procI)
00605 {
00606 sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00607 }
00608 Pstream::gatherList(sendSizes);
00609 Pstream::scatterList(sendSizes);
00610
00611
00612
00613
00614
00615 labelListList constructMap(Pstream::nProcs());
00616
00617
00618 constructMap[Pstream::myProcNo()] = identity
00619 (
00620 sendMap[Pstream::myProcNo()].size()
00621 );
00622
00623 label segmentI = constructMap[Pstream::myProcNo()].size();
00624 forAll(constructMap, procI)
00625 {
00626 if (procI != Pstream::myProcNo())
00627 {
00628
00629 label nRecv = sendSizes[procI][Pstream::myProcNo()];
00630 constructMap[procI].setSize(nRecv);
00631
00632 for (label i = 0; i < nRecv; i++)
00633 {
00634 constructMap[procI][i] = segmentI++;
00635 }
00636 }
00637 }
00638
00639
00640
00641
00642
00643 autoPtr<mapDistribute> mapPtr
00644 (
00645 new mapDistribute
00646 (
00647 segmentI,
00648 sendMap,
00649 constructMap,
00650 true
00651 )
00652 );
00653 const mapDistribute& map = mapPtr();
00654
00655
00656
00657
00658
00659 map.distribute
00660 (
00661
00662
00663 Pstream::nonBlocking,
00664 List<labelPair>(0),
00665 map.constructSize(),
00666 map.subMap(),
00667 map.constructMap(),
00668 triangleIndex
00669 );
00670
00671
00672 return mapPtr;
00673 }
00674
00675
00676 Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
00677 (
00678 const point& centre,
00679 const scalar radiusSqr,
00680 boolList& overlaps
00681 ) const
00682 {
00683 overlaps = false;
00684 label nOverlaps = 0;
00685
00686 forAll(procBb_, procI)
00687 {
00688 const List<treeBoundBox>& bbs = procBb_[procI];
00689
00690 forAll(bbs, bbI)
00691 {
00692 if (bbs[bbI].overlaps(centre, radiusSqr))
00693 {
00694 overlaps[procI] = true;
00695 nOverlaps++;
00696 break;
00697 }
00698 }
00699 }
00700 return nOverlaps;
00701 }
00702
00703
00704
00705
00706
00707 Foam::autoPtr<Foam::mapDistribute>
00708 Foam::distributedTriSurfaceMesh::calcLocalQueries
00709 (
00710 const pointField& centres,
00711 const scalarField& radiusSqr,
00712
00713 pointField& allCentres,
00714 scalarField& allRadiusSqr,
00715 labelList& allSegmentMap
00716 ) const
00717 {
00718
00719
00720
00721 labelListList sendMap(Pstream::nProcs());
00722
00723 {
00724
00725 DynamicList<point> dynAllCentres(centres.size());
00726 DynamicList<scalar> dynAllRadiusSqr(centres.size());
00727
00728 DynamicList<label> dynAllSegmentMap(centres.size());
00729
00730 List<DynamicList<label> > dynSendMap(Pstream::nProcs());
00731
00732
00733 boolList procBbOverlaps(Pstream::nProcs());
00734
00735 forAll(centres, centreI)
00736 {
00737
00738 calcOverlappingProcs
00739 (
00740 centres[centreI],
00741 radiusSqr[centreI],
00742 procBbOverlaps
00743 );
00744
00745 forAll(procBbOverlaps, procI)
00746 {
00747 if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
00748 {
00749 dynSendMap[procI].append(dynAllCentres.size());
00750 dynAllSegmentMap.append(centreI);
00751 dynAllCentres.append(centres[centreI]);
00752 dynAllRadiusSqr.append(radiusSqr[centreI]);
00753 }
00754 }
00755 }
00756
00757
00758 sendMap.setSize(Pstream::nProcs());
00759 forAll(sendMap, procI)
00760 {
00761 dynSendMap[procI].shrink();
00762 sendMap[procI].transfer(dynSendMap[procI]);
00763 }
00764
00765 allCentres.transfer(dynAllCentres.shrink());
00766 allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
00767 allSegmentMap.transfer(dynAllSegmentMap.shrink());
00768 }
00769
00770
00771
00772 labelListList sendSizes(Pstream::nProcs());
00773 sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
00774 forAll(sendMap, procI)
00775 {
00776 sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
00777 }
00778 Pstream::gatherList(sendSizes);
00779 Pstream::scatterList(sendSizes);
00780
00781
00782
00783 labelListList constructMap(Pstream::nProcs());
00784
00785
00786 constructMap[Pstream::myProcNo()] = identity
00787 (
00788 sendMap[Pstream::myProcNo()].size()
00789 );
00790
00791 label segmentI = constructMap[Pstream::myProcNo()].size();
00792 forAll(constructMap, procI)
00793 {
00794 if (procI != Pstream::myProcNo())
00795 {
00796
00797 label nRecv = sendSizes[procI][Pstream::myProcNo()];
00798 constructMap[procI].setSize(nRecv);
00799
00800 for (label i = 0; i < nRecv; i++)
00801 {
00802 constructMap[procI][i] = segmentI++;
00803 }
00804 }
00805 }
00806
00807 autoPtr<mapDistribute> mapPtr
00808 (
00809 new mapDistribute
00810 (
00811 segmentI,
00812 sendMap,
00813 constructMap,
00814 true
00815 )
00816 );
00817 return mapPtr;
00818 }
00819
00820
00821
00822
00823
00824
00825
00826
00827 Foam::List<Foam::List<Foam::treeBoundBox> >
00828 Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
00829 (
00830 const triSurface& s
00831 )
00832 {
00833 if (!decomposer_.valid())
00834 {
00835
00836
00837 IOdictionary decomposeDict
00838 (
00839 IOobject
00840 (
00841 "decomposeParDict",
00842 searchableSurface::time().system(),
00843 searchableSurface::time(),
00844 IOobject::MUST_READ,
00845 IOobject::NO_WRITE,
00846 false
00847 )
00848 );
00849 decomposer_ = decompositionMethod::New(decomposeDict);
00850
00851 if (!decomposer_().parallelAware())
00852 {
00853 FatalErrorIn
00854 (
00855 "distributedTriSurfaceMesh::independentlyDistributedBbs"
00856 "(const triSurface&)"
00857 ) << "The decomposition method " << decomposer_().typeName
00858 << " does not decompose in parallel."
00859 << " Please choose one that does." << exit(FatalError);
00860 }
00861 }
00862
00863
00864 pointField triCentres(s.size());
00865 forAll (s, triI)
00866 {
00867 triCentres[triI] = s[triI].centre(s.points());
00868 }
00869
00870
00871 labelList distribution(decomposer_->decompose(triCentres));
00872
00873
00874
00875
00876 List<List<treeBoundBox> > bbs(Pstream::nProcs());
00877 forAll(bbs, procI)
00878 {
00879 bbs[procI].setSize(1);
00880
00881 bbs[procI][0].min() = point( VGREAT, VGREAT, VGREAT);
00882 bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
00883 }
00884
00885 forAll (s, triI)
00886 {
00887 point& bbMin = bbs[distribution[triI]][0].min();
00888 point& bbMax = bbs[distribution[triI]][0].max();
00889
00890 const labelledTri& f = s[triI];
00891 const point& p0 = s.points()[f[0]];
00892 const point& p1 = s.points()[f[1]];
00893 const point& p2 = s.points()[f[2]];
00894
00895 bbMin = min(bbMin, p0);
00896 bbMin = min(bbMin, p1);
00897 bbMin = min(bbMin, p2);
00898
00899 bbMax = max(bbMax, p0);
00900 bbMax = max(bbMax, p1);
00901 bbMax = max(bbMax, p2);
00902 }
00903
00904
00905 forAll(bbs, procI)
00906 {
00907 forAll(bbs[procI], i)
00908 {
00909 reduce(bbs[procI][i].min(), minOp<point>());
00910 reduce(bbs[procI][i].max(), maxOp<point>());
00911 }
00912 }
00913 return bbs;
00914 }
00915
00916
00917
00918 bool Foam::distributedTriSurfaceMesh::overlaps
00919 (
00920 const List<treeBoundBox>& bbs,
00921 const point& p0,
00922 const point& p1,
00923 const point& p2
00924 )
00925 {
00926 forAll(bbs, bbI)
00927 {
00928 const treeBoundBox& bb = bbs[bbI];
00929
00930 treeBoundBox triBb(p0, p0);
00931 triBb.min() = min(triBb.min(), p1);
00932 triBb.min() = min(triBb.min(), p2);
00933
00934 triBb.max() = max(triBb.max(), p1);
00935 triBb.max() = max(triBb.max(), p2);
00936
00937
00938
00939
00940
00941 if (bb.overlaps(triBb))
00942 {
00943
00944 if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
00945 {
00946
00947 return true;
00948 }
00949
00950
00951
00952
00953 bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
00954
00955 if (intersect)
00956 {
00957 return true;
00958 }
00959 }
00960 }
00961 return false;
00962 }
00963
00964
00965 void Foam::distributedTriSurfaceMesh::subsetMeshMap
00966 (
00967 const triSurface& s,
00968 const boolList& include,
00969 const label nIncluded,
00970 labelList& newToOldPoints,
00971 labelList& oldToNewPoints,
00972 labelList& newToOldFaces
00973 )
00974 {
00975 newToOldFaces.setSize(nIncluded);
00976 newToOldPoints.setSize(s.points().size());
00977 oldToNewPoints.setSize(s.points().size());
00978 oldToNewPoints = -1;
00979 {
00980 label faceI = 0;
00981 label pointI = 0;
00982
00983 forAll(include, oldFacei)
00984 {
00985 if (include[oldFacei])
00986 {
00987
00988 newToOldFaces[faceI++] = oldFacei;
00989
00990
00991 const labelledTri& tri = s[oldFacei];
00992
00993 forAll(tri, fp)
00994 {
00995 label oldPointI = tri[fp];
00996
00997 if (oldToNewPoints[oldPointI] == -1)
00998 {
00999 oldToNewPoints[oldPointI] = pointI;
01000 newToOldPoints[pointI++] = oldPointI;
01001 }
01002 }
01003 }
01004 }
01005 newToOldPoints.setSize(pointI);
01006 }
01007 }
01008
01009
01010 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01011 (
01012 const triSurface& s,
01013 const labelList& newToOldPoints,
01014 const labelList& oldToNewPoints,
01015 const labelList& newToOldFaces
01016 )
01017 {
01018
01019 pointField newPoints(newToOldPoints.size());
01020 forAll(newToOldPoints, i)
01021 {
01022 newPoints[i] = s.points()[newToOldPoints[i]];
01023 }
01024
01025 List<labelledTri> newTriangles(newToOldFaces.size());
01026
01027 forAll(newToOldFaces, i)
01028 {
01029
01030 const labelledTri& tri = s[newToOldFaces[i]];
01031
01032 newTriangles[i][0] = oldToNewPoints[tri[0]];
01033 newTriangles[i][1] = oldToNewPoints[tri[1]];
01034 newTriangles[i][2] = oldToNewPoints[tri[2]];
01035 newTriangles[i].region() = tri.region();
01036 }
01037
01038
01039 return triSurface(newTriangles, s.patches(), newPoints, true);
01040 }
01041
01042
01043 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01044 (
01045 const triSurface& s,
01046 const boolList& include,
01047 labelList& newToOldPoints,
01048 labelList& newToOldFaces
01049 )
01050 {
01051 label n = 0;
01052
01053 forAll(include, i)
01054 {
01055 if (include[i])
01056 {
01057 n++;
01058 }
01059 }
01060
01061 labelList oldToNewPoints;
01062 subsetMeshMap
01063 (
01064 s,
01065 include,
01066 n,
01067 newToOldPoints,
01068 oldToNewPoints,
01069 newToOldFaces
01070 );
01071
01072 return subsetMesh
01073 (
01074 s,
01075 newToOldPoints,
01076 oldToNewPoints,
01077 newToOldFaces
01078 );
01079 }
01080
01081
01082 Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
01083 (
01084 const triSurface& s,
01085 const labelList& newToOldFaces,
01086 labelList& newToOldPoints
01087 )
01088 {
01089 const boolList include
01090 (
01091 createWithValues<boolList>
01092 (
01093 s.size(),
01094 false,
01095 newToOldFaces,
01096 true
01097 )
01098 );
01099
01100 newToOldPoints.setSize(s.points().size());
01101 labelList oldToNewPoints(s.points().size(), -1);
01102 {
01103 label pointI = 0;
01104
01105 forAll(include, oldFacei)
01106 {
01107 if (include[oldFacei])
01108 {
01109
01110 const labelledTri& tri = s[oldFacei];
01111
01112 forAll(tri, fp)
01113 {
01114 label oldPointI = tri[fp];
01115
01116 if (oldToNewPoints[oldPointI] == -1)
01117 {
01118 oldToNewPoints[oldPointI] = pointI;
01119 newToOldPoints[pointI++] = oldPointI;
01120 }
01121 }
01122 }
01123 }
01124 newToOldPoints.setSize(pointI);
01125 }
01126
01127 return subsetMesh
01128 (
01129 s,
01130 newToOldPoints,
01131 oldToNewPoints,
01132 newToOldFaces
01133 );
01134 }
01135
01136
01137 Foam::label Foam::distributedTriSurfaceMesh::findTriangle
01138 (
01139 const List<labelledTri>& allFaces,
01140 const labelListList& allPointFaces,
01141 const labelledTri& otherF
01142 )
01143 {
01144
01145 const labelList& pFaces = allPointFaces[otherF[0]];
01146
01147 forAll(pFaces, i)
01148 {
01149 const labelledTri& f = allFaces[pFaces[i]];
01150
01151 if (f.region() == otherF.region())
01152 {
01153
01154 label fp0 = findIndex(f, otherF[0]);
01155
01156 label fp1 = f.fcIndex(fp0);
01157 label fp2 = f.fcIndex(fp1);
01158
01159 if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
01160 {
01161 return pFaces[i];
01162 }
01163 }
01164 }
01165 return -1;
01166 }
01167
01168
01169
01170 void Foam::distributedTriSurfaceMesh::merge
01171 (
01172 const scalar mergeDist,
01173 const List<labelledTri>& subTris,
01174 const pointField& subPoints,
01175
01176 List<labelledTri>& allTris,
01177 pointField& allPoints,
01178
01179 labelList& faceConstructMap,
01180 labelList& pointConstructMap
01181 )
01182 {
01183 labelList subToAll;
01184 matchPoints
01185 (
01186 subPoints,
01187 allPoints,
01188 scalarField(subPoints.size(), mergeDist),
01189 false,
01190 pointConstructMap
01191 );
01192
01193 label nOldAllPoints = allPoints.size();
01194
01195
01196
01197
01198
01199 label allPointI = nOldAllPoints;
01200 forAll(pointConstructMap, pointI)
01201 {
01202 if (pointConstructMap[pointI] == -1)
01203 {
01204 pointConstructMap[pointI] = allPointI++;
01205 }
01206 }
01207
01208 if (allPointI > nOldAllPoints)
01209 {
01210 allPoints.setSize(allPointI);
01211
01212 forAll(pointConstructMap, pointI)
01213 {
01214 if (pointConstructMap[pointI] >= nOldAllPoints)
01215 {
01216 allPoints[pointConstructMap[pointI]] = subPoints[pointI];
01217 }
01218 }
01219 }
01220
01221
01222
01223 labelListList allPointFaces;
01224 invertManyToMany(nOldAllPoints, allTris, allPointFaces);
01225
01226
01227
01228
01229
01230 label allTriI = allTris.size();
01231 allTris.setSize(allTriI + subTris.size());
01232
01233 faceConstructMap.setSize(subTris.size());
01234
01235 forAll(subTris, triI)
01236 {
01237 const labelledTri& subTri = subTris[triI];
01238
01239
01240 labelledTri mappedTri
01241 (
01242 pointConstructMap[subTri[0]],
01243 pointConstructMap[subTri[1]],
01244 pointConstructMap[subTri[2]],
01245 subTri.region()
01246 );
01247
01248
01249
01250 bool fullMatch = true;
01251
01252 forAll(mappedTri, fp)
01253 {
01254 if (mappedTri[fp] >= nOldAllPoints)
01255 {
01256 fullMatch = false;
01257 break;
01258 }
01259 }
01260
01261 if (fullMatch)
01262 {
01263
01264
01265 label i = findTriangle
01266 (
01267 allTris,
01268 allPointFaces,
01269 mappedTri
01270 );
01271
01272 if (i == -1)
01273 {
01274
01275 faceConstructMap[triI] = allTriI;
01276 allTris[allTriI] = mappedTri;
01277 allTriI++;
01278 }
01279 else
01280 {
01281 faceConstructMap[triI] = i;
01282 }
01283 }
01284 else
01285 {
01286
01287 faceConstructMap[triI] = allTriI;
01288 allTris[allTriI] = mappedTri;
01289 allTriI++;
01290 }
01291 }
01292 allTris.setSize(allTriI);
01293 }
01294
01295
01296
01297
01298 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
01299 (
01300 const IOobject& io,
01301 const triSurface& s,
01302 const dictionary& dict
01303 )
01304 :
01305 triSurfaceMesh(io, s),
01306 dict_
01307 (
01308 IOobject
01309 (
01310 searchableSurface::name() + "Dict",
01311 searchableSurface::instance(),
01312 searchableSurface::local(),
01313 searchableSurface::db(),
01314 searchableSurface::NO_READ,
01315 searchableSurface::writeOpt(),
01316 searchableSurface::registerObject()
01317 ),
01318 dict
01319 )
01320 {
01321 read();
01322
01323 if (debug)
01324 {
01325 Info<< "Constructed from triSurface:" << endl;
01326 writeStats(Info);
01327
01328 labelList nTris(Pstream::nProcs());
01329 nTris[Pstream::myProcNo()] = triSurface::size();
01330 Pstream::gatherList(nTris);
01331 Pstream::scatterList(nTris);
01332
01333 Info<< endl<< "\tproc\ttris\tbb" << endl;
01334 forAll(nTris, procI)
01335 {
01336 Info<< '\t' << procI << '\t' << nTris[procI]
01337 << '\t' << procBb_[procI] << endl;
01338 }
01339 Info<< endl;
01340 }
01341 }
01342
01343
01344 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
01345 :
01346
01347 triSurfaceMesh
01348 (
01349 IOobject
01350 (
01351 io.name(),
01352 io.time().findInstance(io.local(), word::null),
01353 io.local(),
01354 io.db(),
01355 io.readOpt(),
01356 io.writeOpt(),
01357 io.registerObject()
01358 )
01359 ),
01360 dict_
01361 (
01362 IOobject
01363 (
01364 searchableSurface::name() + "Dict",
01365 searchableSurface::instance(),
01366 searchableSurface::local(),
01367 searchableSurface::db(),
01368 searchableSurface::readOpt(),
01369 searchableSurface::writeOpt(),
01370 searchableSurface::registerObject()
01371 )
01372 )
01373 {
01374 read();
01375
01376 if (debug)
01377 {
01378 Info<< "Read distributedTriSurface from " << io.objectPath()
01379 << ':' << endl;
01380 writeStats(Info);
01381
01382 labelList nTris(Pstream::nProcs());
01383 nTris[Pstream::myProcNo()] = triSurface::size();
01384 Pstream::gatherList(nTris);
01385 Pstream::scatterList(nTris);
01386
01387 Info<< endl<< "\tproc\ttris\tbb" << endl;
01388 forAll(nTris, procI)
01389 {
01390 Info<< '\t' << procI << '\t' << nTris[procI]
01391 << '\t' << procBb_[procI] << endl;
01392 }
01393 Info<< endl;
01394 }
01395 }
01396
01397
01398 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
01399 (
01400 const IOobject& io,
01401 const dictionary& dict
01402 )
01403 :
01404
01405 triSurfaceMesh
01406 (
01407 IOobject
01408 (
01409 io.name(),
01410 io.time().findInstance(io.local(), word::null),
01411 io.local(),
01412 io.db(),
01413 io.readOpt(),
01414 io.writeOpt(),
01415 io.registerObject()
01416 ),
01417 dict
01418 ),
01419 dict_
01420 (
01421 IOobject
01422 (
01423 searchableSurface::name() + "Dict",
01424 searchableSurface::instance(),
01425 searchableSurface::local(),
01426 searchableSurface::db(),
01427 searchableSurface::readOpt(),
01428 searchableSurface::writeOpt(),
01429 searchableSurface::registerObject()
01430 )
01431 )
01432 {
01433 read();
01434
01435 if (debug)
01436 {
01437 Info<< "Read distributedTriSurface from " << io.objectPath()
01438 << " and dictionary:" << endl;
01439 writeStats(Info);
01440
01441 labelList nTris(Pstream::nProcs());
01442 nTris[Pstream::myProcNo()] = triSurface::size();
01443 Pstream::gatherList(nTris);
01444 Pstream::scatterList(nTris);
01445
01446 Info<< endl<< "\tproc\ttris\tbb" << endl;
01447 forAll(nTris, procI)
01448 {
01449 Info<< '\t' << procI << '\t' << nTris[procI]
01450 << '\t' << procBb_[procI] << endl;
01451 }
01452 Info<< endl;
01453 }
01454 }
01455
01456
01457
01458
01459 Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
01460 {
01461 clearOut();
01462 }
01463
01464
01465 void Foam::distributedTriSurfaceMesh::clearOut()
01466 {
01467 globalTris_.clear();
01468 triSurfaceMesh::clearOut();
01469 }
01470
01471
01472
01473
01474 const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
01475 {
01476 if (!globalTris_.valid())
01477 {
01478 globalTris_.reset(new globalIndex(triSurface::size()));
01479 }
01480 return globalTris_;
01481 }
01482
01483
01484 void Foam::distributedTriSurfaceMesh::findNearest
01485 (
01486 const pointField& samples,
01487 const scalarField& nearestDistSqr,
01488 List<pointIndexHit>& info
01489 ) const
01490 {
01491 const indexedOctree<treeDataTriSurface>& octree = tree();
01492
01493
01494 const globalIndex& triIndexer = globalTris();
01495
01496
01497
01498
01499
01500 info.setSize(samples.size());
01501 forAll(info, i)
01502 {
01503 info[i].setMiss();
01504 }
01505
01506
01507
01508
01509
01510
01511 label nLocal = 0;
01512
01513 {
01514
01515 boolList procBbOverlaps(Pstream::nProcs());
01516
01517 forAll(samples, i)
01518 {
01519
01520 label nProcs = calcOverlappingProcs
01521 (
01522 samples[i],
01523 nearestDistSqr[i],
01524 procBbOverlaps
01525 );
01526
01527
01528 if (procBbOverlaps[Pstream::myProcNo()])
01529 {
01530 info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
01531 if (info[i].hit())
01532 {
01533 info[i].setIndex(triIndexer.toGlobal(info[i].index()));
01534 }
01535 if (nProcs == 1)
01536 {
01537
01538 nLocal++;
01539 }
01540 }
01541 }
01542 }
01543
01544
01545 if
01546 (
01547 Pstream::parRun()
01548 && (
01549 returnReduce(nLocal, sumOp<label>())
01550 < returnReduce(samples.size(), sumOp<label>())
01551 )
01552 )
01553 {
01554
01555
01556
01557
01558
01559
01560 pointField allCentres;
01561 scalarField allRadiusSqr;
01562 labelList allSegmentMap;
01563 autoPtr<mapDistribute> mapPtr
01564 (
01565 calcLocalQueries
01566 (
01567 samples,
01568 nearestDistSqr,
01569
01570 allCentres,
01571 allRadiusSqr,
01572 allSegmentMap
01573 )
01574 );
01575 const mapDistribute& map = mapPtr();
01576
01577
01578
01579
01580
01581 map.distribute
01582 (
01583
01584
01585 Pstream::nonBlocking,
01586 List<labelPair>(0),
01587 map.constructSize(),
01588 map.subMap(),
01589 map.constructMap(),
01590 allCentres
01591 );
01592 map.distribute
01593 (
01594
01595
01596 Pstream::nonBlocking,
01597 List<labelPair>(0),
01598 map.constructSize(),
01599 map.subMap(),
01600 map.constructMap(),
01601 allRadiusSqr
01602 );
01603
01604
01605
01606
01607
01608 List<pointIndexHit> allInfo(allCentres.size());
01609 forAll(allInfo, i)
01610 {
01611 allInfo[i] = octree.findNearest
01612 (
01613 allCentres[i],
01614 allRadiusSqr[i]
01615 );
01616 if (allInfo[i].hit())
01617 {
01618 allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
01619 }
01620 }
01621
01622
01623
01624
01625
01626 map.distribute
01627 (
01628
01629
01630
01631
01632
01633
01634 Pstream::nonBlocking,
01635 List<labelPair>(0),
01636 allSegmentMap.size(),
01637 map.constructMap(),
01638 map.subMap(),
01639 allInfo
01640 );
01641
01642
01643
01644
01645
01646 forAll(allInfo, i)
01647 {
01648 if (allInfo[i].hit())
01649 {
01650 label pointI = allSegmentMap[i];
01651
01652 if (!info[pointI].hit())
01653 {
01654
01655 info[pointI] = allInfo[i];
01656 }
01657 else
01658 {
01659
01660 if
01661 (
01662 magSqr(allInfo[i].hitPoint()-samples[pointI])
01663 < magSqr(info[pointI].hitPoint()-samples[pointI])
01664 )
01665 {
01666 info[pointI] = allInfo[i];
01667 }
01668 }
01669 }
01670 }
01671 }
01672 }
01673
01674
01675 void Foam::distributedTriSurfaceMesh::findLine
01676 (
01677 const pointField& start,
01678 const pointField& end,
01679 List<pointIndexHit>& info
01680 ) const
01681 {
01682 findLine
01683 (
01684 true,
01685 start,
01686 end,
01687 info
01688 );
01689 }
01690
01691
01692 void Foam::distributedTriSurfaceMesh::findLineAny
01693 (
01694 const pointField& start,
01695 const pointField& end,
01696 List<pointIndexHit>& info
01697 ) const
01698 {
01699 findLine
01700 (
01701 true,
01702 start,
01703 end,
01704 info
01705 );
01706 }
01707
01708
01709 void Foam::distributedTriSurfaceMesh::findLineAll
01710 (
01711 const pointField& start,
01712 const pointField& end,
01713 List<List<pointIndexHit> >& info
01714 ) const
01715 {
01716
01717
01718
01719
01720
01721
01722 List<pointIndexHit> hitInfo(start.size());
01723
01724 findLine
01725 (
01726 true,
01727 start,
01728 end,
01729 hitInfo
01730 );
01731
01732
01733
01734
01735
01736
01737
01738 const vectorField dirVec(end-start);
01739 const scalarField magSqrDirVec(magSqr(dirVec));
01740 const vectorField smallVec
01741 (
01742 Foam::sqrt(SMALL)*dirVec
01743 + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
01744 );
01745
01746
01747 labelList pointMap(start.size());
01748 pointField e0(start.size());
01749 pointField e1(start.size());
01750 label compactI = 0;
01751
01752 info.setSize(hitInfo.size());
01753 forAll(hitInfo, pointI)
01754 {
01755 if (hitInfo[pointI].hit())
01756 {
01757 info[pointI].setSize(1);
01758 info[pointI][0] = hitInfo[pointI];
01759
01760 point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];
01761
01762 if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
01763 {
01764 e0[compactI] = pt;
01765 e1[compactI] = end[pointI];
01766 pointMap[compactI] = pointI;
01767 compactI++;
01768 }
01769 }
01770 else
01771 {
01772 info[pointI].clear();
01773 }
01774 }
01775
01776 e0.setSize(compactI);
01777 e1.setSize(compactI);
01778 pointMap.setSize(compactI);
01779
01780 while (returnReduce(e0.size(), sumOp<label>()) > 0)
01781 {
01782 findLine
01783 (
01784 true,
01785 e0,
01786 e1,
01787 hitInfo
01788 );
01789
01790
01791
01792 label compactI = 0;
01793 forAll(hitInfo, i)
01794 {
01795 if (hitInfo[i].hit())
01796 {
01797 label pointI = pointMap[i];
01798
01799 label sz = info[pointI].size();
01800 info[pointI].setSize(sz+1);
01801 info[pointI][sz] = hitInfo[i];
01802
01803 point pt = hitInfo[i].hitPoint() + smallVec[pointI];
01804
01805 if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
01806 {
01807 e0[compactI] = pt;
01808 e1[compactI] = end[pointI];
01809 pointMap[compactI] = pointI;
01810 compactI++;
01811 }
01812 }
01813 }
01814
01815
01816 e0.setSize(compactI);
01817 e1.setSize(compactI);
01818 pointMap.setSize(compactI);
01819 }
01820 }
01821
01822
01823 void Foam::distributedTriSurfaceMesh::getRegion
01824 (
01825 const List<pointIndexHit>& info,
01826 labelList& region
01827 ) const
01828 {
01829 if (!Pstream::parRun())
01830 {
01831 region.setSize(info.size());
01832 forAll(info, i)
01833 {
01834 if (info[i].hit())
01835 {
01836 region[i] = triSurface::operator[](info[i].index()).region();
01837 }
01838 else
01839 {
01840 region[i] = -1;
01841 }
01842 }
01843 return;
01844 }
01845
01846
01847
01848
01849
01850 labelList triangleIndex(info.size());
01851 autoPtr<mapDistribute> mapPtr
01852 (
01853 calcLocalQueries
01854 (
01855 info,
01856 triangleIndex
01857 )
01858 );
01859 const mapDistribute& map = mapPtr();
01860
01861
01862
01863
01864
01865 const triSurface& s = static_cast<const triSurface&>(*this);
01866
01867 region.setSize(triangleIndex.size());
01868
01869 forAll(triangleIndex, i)
01870 {
01871 label triI = triangleIndex[i];
01872 region[i] = s[triI].region();
01873 }
01874
01875
01876
01877
01878
01879 map.distribute
01880 (
01881
01882
01883
01884
01885
01886
01887 Pstream::nonBlocking,
01888 List<labelPair>(0),
01889 info.size(),
01890 map.constructMap(),
01891 map.subMap(),
01892 region
01893 );
01894 }
01895
01896
01897 void Foam::distributedTriSurfaceMesh::getNormal
01898 (
01899 const List<pointIndexHit>& info,
01900 vectorField& normal
01901 ) const
01902 {
01903 if (!Pstream::parRun())
01904 {
01905 triSurfaceMesh::getNormal(info, normal);
01906 return;
01907 }
01908
01909
01910
01911
01912
01913 labelList triangleIndex(info.size());
01914 autoPtr<mapDistribute> mapPtr
01915 (
01916 calcLocalQueries
01917 (
01918 info,
01919 triangleIndex
01920 )
01921 );
01922 const mapDistribute& map = mapPtr();
01923
01924
01925
01926
01927
01928 const triSurface& s = static_cast<const triSurface&>(*this);
01929
01930 normal.setSize(triangleIndex.size());
01931
01932 forAll(triangleIndex, i)
01933 {
01934 label triI = triangleIndex[i];
01935 normal[i] = s[triI].normal(s.points());
01936 normal[i] /= mag(normal[i]) + VSMALL;
01937 }
01938
01939
01940
01941
01942
01943 map.distribute
01944 (
01945
01946
01947
01948
01949
01950
01951 Pstream::nonBlocking,
01952 List<labelPair>(0),
01953 info.size(),
01954 map.constructMap(),
01955 map.subMap(),
01956 normal
01957 );
01958 }
01959
01960
01961 void Foam::distributedTriSurfaceMesh::getField
01962 (
01963 const List<pointIndexHit>& info,
01964 labelList& values
01965 ) const
01966 {
01967 if (!Pstream::parRun())
01968 {
01969 triSurfaceMesh::getField(info, values);
01970 return;
01971 }
01972
01973 if (foundObject<triSurfaceLabelField>("values"))
01974 {
01975 const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
01976 (
01977 "values"
01978 );
01979
01980
01981
01982
01983
01984 labelList triangleIndex(info.size());
01985 autoPtr<mapDistribute> mapPtr
01986 (
01987 calcLocalQueries
01988 (
01989 info,
01990 triangleIndex
01991 )
01992 );
01993 const mapDistribute& map = mapPtr();
01994
01995
01996
01997
01998
01999 values.setSize(triangleIndex.size());
02000
02001 forAll(triangleIndex, i)
02002 {
02003 label triI = triangleIndex[i];
02004 values[i] = fld[triI];
02005 }
02006
02007
02008
02009
02010
02011 map.distribute
02012 (
02013 Pstream::nonBlocking,
02014 List<labelPair>(0),
02015 info.size(),
02016 map.constructMap(),
02017 map.subMap(),
02018 values
02019 );
02020 }
02021 }
02022
02023
02024 void Foam::distributedTriSurfaceMesh::getVolumeType
02025 (
02026 const pointField& points,
02027 List<volumeType>& volType
02028 ) const
02029 {
02030 FatalErrorIn
02031 (
02032 "distributedTriSurfaceMesh::getVolumeType"
02033 "(const pointField&, List<volumeType>&) const"
02034 ) << "Volume type not supported for distributed surfaces."
02035 << exit(FatalError);
02036 }
02037
02038
02039
02040 Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
02041 (
02042 const triSurface& s,
02043 const List<treeBoundBox>& bbs,
02044
02045 labelList& subPointMap,
02046 labelList& subFaceMap
02047 )
02048 {
02049
02050 boolList includedFace(s.size(), false);
02051
02052
02053 List<treeBoundBox> bbsX(bbs.size());
02054 const scalar eps = 1.0e-4;
02055 forAll(bbs, i)
02056 {
02057 const point mid = 0.5*(bbs[i].min() + bbs[i].max());
02058 const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
02059
02060 bbsX[i].min() = mid - halfSpan;
02061 bbsX[i].max() = mid + halfSpan;
02062 }
02063
02064 forAll(s, triI)
02065 {
02066 const labelledTri& f = s[triI];
02067 const point& p0 = s.points()[f[0]];
02068 const point& p1 = s.points()[f[1]];
02069 const point& p2 = s.points()[f[2]];
02070
02071 if (overlaps(bbsX, p0, p1, p2))
02072 {
02073 includedFace[triI] = true;
02074 }
02075 }
02076
02077 return subsetMesh(s, includedFace, subPointMap, subFaceMap);
02078 }
02079
02080
02081 void Foam::distributedTriSurfaceMesh::distribute
02082 (
02083 const List<treeBoundBox>& bbs,
02084 const bool keepNonLocal,
02085 autoPtr<mapDistribute>& faceMap,
02086 autoPtr<mapDistribute>& pointMap
02087 )
02088 {
02089
02090
02091
02092 {
02093 List<List<treeBoundBox> > newProcBb(Pstream::nProcs());
02094
02095 switch(distType_)
02096 {
02097 case FOLLOW:
02098 newProcBb[Pstream::myProcNo()].setSize(bbs.size());
02099 forAll(bbs, i)
02100 {
02101 newProcBb[Pstream::myProcNo()][i] = bbs[i];
02102 }
02103 Pstream::gatherList(newProcBb);
02104 Pstream::scatterList(newProcBb);
02105 break;
02106
02107 case INDEPENDENT:
02108 newProcBb = independentlyDistributedBbs(*this);
02109 break;
02110
02111 case FROZEN:
02112 return;
02113 break;
02114
02115 default:
02116 FatalErrorIn("distributedTriSurfaceMesh::distribute(..)")
02117 << "Unsupported distribution type." << exit(FatalError);
02118 break;
02119 }
02120
02121
02122
02123
02124
02125
02126
02127
02128 if (newProcBb == procBb_)
02129 {
02130 return;
02131 }
02132 else
02133 {
02134 procBb_.transfer(newProcBb);
02135 dict_.set("bounds", procBb_[Pstream::myProcNo()]);
02136 }
02137 }
02138
02139
02140
02141 if (debug)
02142 {
02143 labelList nTris(Pstream::nProcs());
02144 nTris[Pstream::myProcNo()] = triSurface::size();
02145 Pstream::gatherList(nTris);
02146 Pstream::scatterList(nTris);
02147
02148 Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
02149 << endl
02150 << "\tproc\ttris" << endl;
02151
02152 forAll(nTris, procI)
02153 {
02154 Info<< '\t' << procI << '\t' << nTris[procI] << endl;
02155 }
02156 Info<< endl;
02157 }
02158
02159
02160
02161
02162
02163 labelListList faceSendMap(Pstream::nProcs());
02164 labelListList pointSendMap(Pstream::nProcs());
02165
02166 forAll(procBb_, procI)
02167 {
02168 overlappingSurface
02169 (
02170 *this,
02171 procBb_[procI],
02172 pointSendMap[procI],
02173 faceSendMap[procI]
02174 );
02175
02176 if (debug)
02177 {
02178
02179
02180
02181 }
02182 }
02183
02184 if (keepNonLocal)
02185 {
02186
02187
02188
02189 const triSurface& s = static_cast<const triSurface&>(*this);
02190
02191 boolList includedFace(s.size(), true);
02192
02193 forAll(faceSendMap, procI)
02194 {
02195 if (procI != Pstream::myProcNo())
02196 {
02197 forAll(faceSendMap[procI], i)
02198 {
02199 includedFace[faceSendMap[procI][i]] = false;
02200 }
02201 }
02202 }
02203
02204
02205
02206
02207 forAll(faceSendMap[Pstream::myProcNo()], i)
02208 {
02209 includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
02210 }
02211
02212 subsetMesh
02213 (
02214 s,
02215 includedFace,
02216 pointSendMap[Pstream::myProcNo()],
02217 faceSendMap[Pstream::myProcNo()]
02218 );
02219 }
02220
02221
02222
02223
02224
02225 labelListList faceSendSizes(Pstream::nProcs());
02226 faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
02227 forAll(faceSendMap, procI)
02228 {
02229 faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
02230 }
02231 Pstream::gatherList(faceSendSizes);
02232 Pstream::scatterList(faceSendSizes);
02233
02234
02235
02236
02237
02238
02239
02240 List<labelledTri> allTris;
02241 pointField allPoints;
02242
02243 labelListList faceConstructMap(Pstream::nProcs());
02244 labelListList pointConstructMap(Pstream::nProcs());
02245
02246
02247
02248
02249
02250 {
02251 labelList pointMap;
02252 triSurface subSurface
02253 (
02254 subsetMesh
02255 (
02256 *this,
02257 faceSendMap[Pstream::myProcNo()],
02258 pointMap
02259 )
02260 );
02261
02262 allTris = subSurface;
02263 allPoints = subSurface.points();
02264
02265 faceConstructMap[Pstream::myProcNo()] = identity
02266 (
02267 faceSendMap[Pstream::myProcNo()].size()
02268 );
02269 pointConstructMap[Pstream::myProcNo()] = identity
02270 (
02271 pointSendMap[Pstream::myProcNo()].size()
02272 );
02273 }
02274
02275
02276
02277
02278
02279
02280 forAll(faceSendSizes, procI)
02281 {
02282 if (procI != Pstream::myProcNo())
02283 {
02284 if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
02285 {
02286 OPstream str(Pstream::blocking, procI);
02287
02288 labelList pointMap;
02289 triSurface subSurface
02290 (
02291 subsetMesh
02292 (
02293 *this,
02294 faceSendMap[procI],
02295 pointMap
02296 )
02297 );
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307 str << subSurface;
02308 }
02309 }
02310 }
02311
02312
02313
02314
02315
02316 forAll(faceSendSizes, procI)
02317 {
02318 if (procI != Pstream::myProcNo())
02319 {
02320 if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
02321 {
02322 IPstream str(Pstream::blocking, procI);
02323
02324
02325 triSurface subSurface(str);
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336 merge
02337 (
02338 mergeDist_,
02339 subSurface,
02340 subSurface.points(),
02341
02342 allTris,
02343 allPoints,
02344 faceConstructMap[procI],
02345 pointConstructMap[procI]
02346 );
02347
02348
02349
02350
02351
02352
02353 }
02354 }
02355 }
02356
02357
02358 faceMap.reset
02359 (
02360 new mapDistribute
02361 (
02362 allTris.size(),
02363 faceSendMap,
02364 faceConstructMap,
02365 true
02366 )
02367 );
02368 pointMap.reset
02369 (
02370 new mapDistribute
02371 (
02372 allPoints.size(),
02373 pointSendMap,
02374 pointConstructMap,
02375 true
02376 )
02377 );
02378
02379
02380 triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
02381
02382 clearOut();
02383
02384
02385
02386
02387 distributeFields<label>(faceMap());
02388 distributeFields<scalar>(faceMap());
02389 distributeFields<vector>(faceMap());
02390 distributeFields<sphericalTensor>(faceMap());
02391 distributeFields<symmTensor>(faceMap());
02392 distributeFields<tensor>(faceMap());
02393
02394 if (debug)
02395 {
02396 labelList nTris(Pstream::nProcs());
02397 nTris[Pstream::myProcNo()] = triSurface::size();
02398 Pstream::gatherList(nTris);
02399 Pstream::scatterList(nTris);
02400
02401 Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
02402 << endl
02403 << "\tproc\ttris" << endl;
02404
02405 forAll(nTris, procI)
02406 {
02407 Info<< '\t' << procI << '\t' << nTris[procI] << endl;
02408 }
02409 Info<< endl;
02410 }
02411 }
02412
02413
02414
02415 bool Foam::distributedTriSurfaceMesh::writeObject
02416 (
02417 IOstream::streamFormat fmt,
02418 IOstream::versionNumber ver,
02419 IOstream::compressionType cmp
02420 ) const
02421 {
02422
02423 const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
02424
02425
02426 return
02427 triSurfaceMesh::writeObject(fmt, ver, cmp)
02428 && dict_.writeObject(IOstream::ASCII, ver, cmp);
02429 }
02430
02431
02432 void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
02433 {
02434 boundBox bb;
02435 label nPoints;
02436 calcBounds(bb, nPoints);
02437 reduce(bb.min(), minOp<point>());
02438 reduce(bb.max(), maxOp<point>());
02439
02440 os << "Triangles : " << returnReduce(triSurface::size(), sumOp<label>())
02441 << endl
02442 << "Vertices : " << returnReduce(nPoints, sumOp<label>()) << endl
02443 << "Bounding Box : " << bb << endl;
02444 }
02445
02446
02447