shootRays.H
Go to the documentation of this file.
1 // All rays expressed as start face (local) index end end face (global)
2 // Pre-size by assuming a certain percentage is visible.
3 
4 // Maximum length for dynamicList
6 (
7  viewFactorDict.lookupOrDefault<label>("maxDynListLength", 100000)
8 );
9 
10 for (label proci = 0; proci < Pstream::nProcs(); proci++)
11 {
12  // Shoot rays from me to proci. Note that even if processor has
13  // 0 faces we still need to call findLine to keep calls synced.
14 
15  DynamicField<point> start(coarseMesh.nFaces());
16  DynamicField<point> end(start.size());
17  DynamicList<label> startIndex(start.size());
18  DynamicList<label> endIndex(start.size());
19 
20  DynamicList<label> startAgg(start.size());
21  DynamicList<label> endAgg(start.size());
22 
23 
24  const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
25  const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
26  const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
27 
28  const pointField& remoteArea = remoteCoarseSf[proci];
29  const pointField& remoteFc = remoteCoarseCf[proci];
30  const labelField& remoteAgg = remoteCoarseAgg[proci];
31 
32  label i = 0;
33  label j = 0;
34  do
35  {
36  for (; i < myFc.size(); i++)
37  {
38  const point& fc = myFc[i];
39  const vector& fA = myArea[i];
40  const label& fAgg = myAgg[i];
41 
42  for (; j < remoteFc.size(); j++)
43  {
44  if (proci != Pstream::myProcNo() || i != j)
45  {
46  const point& remFc = remoteFc[j];
47  const vector& remA = remoteArea[j];
48  const label& remAgg = remoteAgg[j];
49 
50  const vector d(remFc - fc);
51 
52  if (((d & fA) < 0.) && ((d & remA) > 0))
53  {
54  start.append(fc + 0.001*d);
55  startIndex.append(i);
56  startAgg.append(globalNumbering.toGlobal(proci, fAgg));
57  end.append(fc + 0.999*d);
58  label globalI = globalNumbering.toGlobal(proci, j);
59  endIndex.append(globalI);
60  endAgg.append(globalNumbering.toGlobal(proci, remAgg));
61  if (startIndex.size() > maxDynListLength)
62  {
64  << "Dynamic list need from capacity."
65  << "Actual size maxDynListLength : "
67  << abort(FatalError);
68  }
69  }
70  }
71  }
72 
73  if (j == remoteFc.size())
74  {
75  j = 0;
76  }
77  }
78 
79  } while (returnReduce(i < myFc.size(), orOp<bool>()));
80 
81  List<pointIndexHit> hitInfo(startIndex.size());
82  surfacesMesh.findLine(start, end, hitInfo);
83 
84  // Return hit global agglo index
85  labelList aggHitIndex;
86  surfacesMesh.getField(hitInfo, aggHitIndex);
87 
88  DynamicList<label> dRayIs;
89 
90  // Collect the rays which don't have obstacle in between rayStartFace
91  // and rayEndFace. If the ray hit itself get stored in dRayIs
92  forAll(hitInfo, rayI)
93  {
94  if (!hitInfo[rayI].hit())
95  {
96  rayStartFace.append(startIndex[rayI]);
97  rayEndFace.append(endIndex[rayI]);
98  }
99  else if (aggHitIndex[rayI] == startAgg[rayI])
100  {
101  dRayIs.append(rayI);
102  }
103  }
104 
105  start.clear();
106 
107 
108  // Continue rays which hit themself. If they hit the target
109  // agglomeration are added to rayStartFace and rayEndFace
110 
111  bool firstLoop = true;
112  DynamicField<point> startHitItself;
113  DynamicField<point> endHitItself;
114  label iter = 0;
115 
116  do
117  {
118  labelField rayIs;
119  rayIs.transfer(dRayIs);
120  dRayIs.clear();
121  forAll(rayIs, rayI)
122  {
123  const label rayID = rayIs[rayI];
124  label hitIndex = -1;
125 
126  if (firstLoop)
127  {
128  hitIndex = rayIs[rayI];
129  }
130  else
131  {
132  hitIndex = rayI;
133  }
134 
135  if (hitInfo[hitIndex].hit())
136  {
137  if (aggHitIndex[hitIndex] == startAgg[rayID])
138  {
139  const vector& endP = end[rayID];
140  const vector& startP = hitInfo[hitIndex].hitPoint();
141  const vector d(endP - startP);
142 
143  startHitItself.append(startP + 0.01*d);
144  endHitItself.append(startP + 1.01*d);
145 
146  dRayIs.append(rayID);
147  }
148  else if (aggHitIndex[hitIndex] == endAgg[rayID])
149  {
150  rayStartFace.append(startIndex[rayID]);
151  rayEndFace.append(endIndex[rayID]);
152  }
153 
154  }
155  }
156 
157  hitInfo.clear();
158  hitInfo.resize(dRayIs.size());
159 
160  surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
161 
162  surfacesMesh.getField(hitInfo, aggHitIndex);
163 
164 
165  endHitItself.clear();
166  startHitItself.clear();
167  firstLoop = false;
168  iter ++;
169 
170  } while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
171 
172  startIndex.clear();
173  end.clear();
174  endIndex.clear();
175  startAgg.clear();
176  endAgg.clear();
177 }
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
volVectorField vectorField(fieldObject, mesh)
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManip< error > abort(error &err)
Definition: errorManip.H:131
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
error FatalError
distributedTriSurfaceMesh surfacesMesh(IOobject("wallSurface.stl", runTime.constant(), searchableSurface::geometryDir(runTime), runTime, IOobject::NO_READ, IOobject::NO_WRITE), localSurface, dict)
const label maxDynListLength(viewFactorDict.lookupOrDefault< label >("maxDynListLength", 100000))