7for (
const int proci : Pstream::allProcs())
9 std::vector<Point> start;
10 start.reserve(coarseMesh.nFaces());
12 std::vector<Point> end(start.size());
13 end.reserve(start.size());
15 DynamicList<label> startIndex(start.size());
16 DynamicList<label> endIndex(start.size());
18 DynamicList<label> startAgg(start.size());
19 DynamicList<label> endAgg(start.size());
21 const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
22 const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
23 const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
25 const pointField& remoteArea = remoteCoarseSf[proci];
26 const pointField& remoteFc = remoteCoarseCf[proci];
27 const labelField& remoteAgg = remoteCoarseAgg[proci];
33 for (; i < myFc.size(); i++)
35 const point& fc = myFc[i];
36 const vector& fA = myArea[i];
37 const label& fAgg = myAgg[i];
39 for (; j < remoteFc.size(); j++)
41 if (proci != Pstream::myProcNo() || i != j)
43 const point& remFc = remoteFc[j];
44 const vector& remA = remoteArea[j];
45 const label& remAgg = remoteAgg[j];
47 const vector d(remFc - fc);
50 const vector nd = d/mag(d);
51 const vector nfA = fA/mag(fA);
52 const vector nremA = remA/mag(remA);
54 if (((nd & nfA) < 0) && ((nd & nremA) > 0))
56 vector direction(d[0], d[1], d[2]);
58 vector rayEnd(
s + (1-intTol)*direction);
59 end.push_back(
Point(rayEnd[0], rayEnd[1], rayEnd[2]));
61 s +=
vector(intTol*d[0], intTol*d[1], intTol*d[2]);
63 start.push_back(
Point(
s[0],
s[1],
s[2]));
69 globalNumbering.toGlobal(proci, fAgg)
73 globalNumbering.toGlobal(proci, remAgg)
77 label globalI = globalNumbering.toGlobal(proci, j);
78 endIndex.append(globalI);
83 <<
"Dynamic list need from capacity."
84 <<
"Actual size maxDynListLength : "
92 if (j == remoteFc.size())
98 }
while (i < myFc.size());
100 label totalRays(startIndex.size());
101 reduce(totalRays, sumOp<scalar>());
105 Pout<<
"Number of rays :" << totalRays << endl;
108 for (
unsigned long rayI = 0; rayI < start.size(); ++rayI)
110 Segment ray(start[rayI], end[rayI]);
111 Segment_intersection intersects =
tree.any_intersection(ray);
115 rayStartFace.append(startIndex[rayI]);
116 rayEndFace.append(endIndex[rayI]);
123 Pout <<
"hits : " << rayStartFace.size()<< endl;
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))