AMIInterpolationParallelOps.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "AMIInterpolation.H"
27 #include "mergePoints.H"
28 #include "mapDistribute.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 Foam::label Foam::AMIInterpolation::calcDistribution
33 (
34  const primitivePatch& srcPatch,
35  const primitivePatch& tgtPatch
36 ) const
37 {
38  label proci = 0;
39 
40  if (Pstream::parRun())
41  {
42  List<label> facesPresentOnProc(Pstream::nProcs(), 0);
43  if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
44  {
45  facesPresentOnProc[Pstream::myProcNo()] = 1;
46  }
47  else
48  {
49  facesPresentOnProc[Pstream::myProcNo()] = 0;
50  }
51 
52  Pstream::gatherList(facesPresentOnProc);
53  Pstream::scatterList(facesPresentOnProc);
54 
55  label nHaveFaces = sum(facesPresentOnProc);
56 
57  if (nHaveFaces > 1)
58  {
59  proci = -1;
60  if (debug)
61  {
63  << "AMI split across multiple processors" << endl;
64  }
65  }
66  else if (nHaveFaces == 1)
67  {
68  proci = findIndex(facesPresentOnProc, 1);
69  if (debug)
70  {
72  << "AMI local to processor" << proci << endl;
73  }
74  }
75  }
76 
77 
78  // Either not parallel or no faces on any processor
79  return proci;
80 }
81 
82 
83 Foam::label Foam::AMIInterpolation::calcOverlappingProcs
84 (
85  const List<treeBoundBoxList>& procBb,
86  const treeBoundBox& bb,
87  boolList& overlaps
88 ) const
89 {
90  overlaps.setSize(procBb.size());
91  overlaps = false;
92 
93  label nOverlaps = 0;
94 
95  forAll(procBb, proci)
96  {
97  const treeBoundBoxList& bbs = procBb[proci];
98 
99  forAll(bbs, bbI)
100  {
101  if (bbs[bbI].overlaps(bb))
102  {
103  overlaps[proci] = true;
104  nOverlaps++;
105  break;
106  }
107  }
108  }
109  return nOverlaps;
110 }
111 
112 
113 void Foam::AMIInterpolation::distributePatches
114 (
115  const mapDistribute& map,
116  const primitivePatch& pp,
117  const globalIndex& gi,
118  List<faceList>& faces,
119  List<pointField>& points,
120  List<labelList>& faceIDs
121 ) const
122 {
123  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
124 
125  for (label domain = 0; domain < Pstream::nProcs(); domain++)
126  {
127  const labelList& sendElems = map.subMap()[domain];
128 
129  if (domain != Pstream::myProcNo() && sendElems.size())
130  {
131  labelList globalElems(sendElems.size());
132  forAll(sendElems, i)
133  {
134  globalElems[i] = gi.toGlobal(sendElems[i]);
135  }
136 
137  faceList subFaces(UIndirectList<face>(pp, sendElems));
138  primitivePatch subPatch
139  (
140  SubList<face>(subFaces, subFaces.size()),
141  pp.points()
142  );
143 
144  if (debug & 2)
145  {
146  Pout<< "distributePatches: to processor " << domain
147  << " sending faces " << subPatch.faceCentres() << endl;
148  }
149 
150  UOPstream toDomain(domain, pBufs);
151  toDomain
152  << subPatch.localFaces() << subPatch.localPoints()
153  << globalElems;
154  }
155  }
156 
157  // Start receiving
158  pBufs.finishedSends();
159 
160  faces.setSize(Pstream::nProcs());
161  points.setSize(Pstream::nProcs());
162  faceIDs.setSize(Pstream::nProcs());
163 
164  {
165  // Set up 'send' to myself
166  const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
167  faceList subFaces(UIndirectList<face>(pp, sendElems));
168  primitivePatch subPatch
169  (
170  SubList<face>(subFaces, subFaces.size()),
171  pp.points()
172  );
173 
174  // Receive
175  if (debug & 2)
176  {
177  Pout<< "distributePatches: to processor " << Pstream::myProcNo()
178  << " sending faces " << subPatch.faceCentres() << endl;
179  }
180 
181  faces[Pstream::myProcNo()] = subPatch.localFaces();
182  points[Pstream::myProcNo()] = subPatch.localPoints();
183 
184  faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
185  forAll(sendElems, i)
186  {
187  faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
188  }
189  }
190 
191  // Consume
192  for (label domain = 0; domain < Pstream::nProcs(); domain++)
193  {
194  const labelList& recvElems = map.constructMap()[domain];
195 
196  if (domain != Pstream::myProcNo() && recvElems.size())
197  {
198  UIPstream str(domain, pBufs);
199 
200  str >> faces[domain]
201  >> points[domain]
202  >> faceIDs[domain];
203  }
204  }
205 }
206 
207 
208 void Foam::AMIInterpolation::distributeAndMergePatches
209 (
210  const mapDistribute& map,
211  const primitivePatch& tgtPatch,
212  const globalIndex& gi,
213  faceList& tgtFaces,
214  pointField& tgtPoints,
215  labelList& tgtFaceIDs
216 ) const
217 {
218  // Exchange per-processor data
219  List<faceList> allFaces;
220  List<pointField> allPoints;
221  List<labelList> allTgtFaceIDs;
222  distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
223 
224  // Renumber and flatten
225  label nFaces = 0;
226  label nPoints = 0;
227  forAll(allFaces, proci)
228  {
229  nFaces += allFaces[proci].size();
230  nPoints += allPoints[proci].size();
231  }
232 
233  tgtFaces.setSize(nFaces);
234  tgtPoints.setSize(nPoints);
235  tgtFaceIDs.setSize(nFaces);
236 
237  nFaces = 0;
238  nPoints = 0;
239 
240  // My own data first
241  {
242  const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
243  SubList<label>(tgtFaceIDs, faceIDs.size()) = faceIDs;
244 
245  const faceList& fcs = allFaces[Pstream::myProcNo()];
246  forAll(fcs, i)
247  {
248  const face& f = fcs[i];
249  face& newF = tgtFaces[nFaces++];
250  newF.setSize(f.size());
251  forAll(f, fp)
252  {
253  newF[fp] = f[fp] + nPoints;
254  }
255  }
256 
257  const pointField& pts = allPoints[Pstream::myProcNo()];
258  forAll(pts, i)
259  {
260  tgtPoints[nPoints++] = pts[i];
261  }
262  }
263 
264 
265  // Other proc data follows
266  forAll(allFaces, proci)
267  {
268  if (proci != Pstream::myProcNo())
269  {
270  const labelList& faceIDs = allTgtFaceIDs[proci];
271  SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces) = faceIDs;
272 
273  const faceList& fcs = allFaces[proci];
274  forAll(fcs, i)
275  {
276  const face& f = fcs[i];
277  face& newF = tgtFaces[nFaces++];
278  newF.setSize(f.size());
279  forAll(f, fp)
280  {
281  newF[fp] = f[fp] + nPoints;
282  }
283  }
284 
285  const pointField& pts = allPoints[proci];
286  forAll(pts, i)
287  {
288  tgtPoints[nPoints++] = pts[i];
289  }
290  }
291  }
292 
293  // Merge
294  labelList oldToNew;
295  pointField newTgtPoints;
296  bool hasMerged = mergePoints
297  (
298  tgtPoints,
299  small,
300  false,
301  oldToNew,
302  newTgtPoints
303  );
304 
305  if (hasMerged)
306  {
307  if (debug)
308  {
309  Pout<< "Merged from " << tgtPoints.size()
310  << " down to " << newTgtPoints.size() << " points" << endl;
311  }
312 
313  tgtPoints.transfer(newTgtPoints);
314  forAll(tgtFaces, i)
315  {
316  inplaceRenumber(oldToNew, tgtFaces[i]);
317  }
318  }
319 }
320 
321 
323 Foam::AMIInterpolation::calcProcMap
324 (
325  const primitivePatch& srcPatch,
326  const primitivePatch& tgtPatch
327 ) const
328 {
329  // Get decomposition of patch
330  List<treeBoundBoxList> procBb(Pstream::nProcs());
331 
332  if (srcPatch.size())
333  {
335  (
336  1, // For now single bounding box per proc
337  treeBoundBox
338  (
339  srcPatch.points(),
340  srcPatch.meshPoints()
341  )
342  );
343  }
344  else
345  {
346  procBb[Pstream::myProcNo()] = treeBoundBoxList();
347  }
348 
349  // slightly increase size of bounding boxes to allow for cases where
350  // bounding boxes are perfectly aligned
351  forAll(procBb[Pstream::myProcNo()], bbI)
352  {
353  treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
354  bb.inflate(0.01);
355  }
356 
357  Pstream::gatherList(procBb);
358  Pstream::scatterList(procBb);
359 
360 
361  if (debug)
362  {
363  Info<< "Determining extent of srcPatch per processor:" << nl
364  << "\tproc\tbb" << endl;
365  forAll(procBb, proci)
366  {
367  Info<< '\t' << proci << '\t' << procBb[proci] << endl;
368  }
369  }
370 
371 
372  // Determine which faces of tgtPatch overlaps srcPatch per proc
373  const faceList& faces = tgtPatch.localFaces();
374  const pointField& points = tgtPatch.localPoints();
375 
376  labelListList sendMap;
377 
378  {
379  // Per processor indices into all segments to send
380  List<DynamicList<label>> dynSendMap(Pstream::nProcs());
381 
382  // Work array - whether processor bb overlaps the face bounds
383  boolList procBbOverlaps(Pstream::nProcs());
384 
385  forAll(faces, facei)
386  {
387  if (faces[facei].size())
388  {
389  treeBoundBox faceBb(points, faces[facei]);
390 
391  // Find the processor this face overlaps
392  calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
393 
394  forAll(procBbOverlaps, proci)
395  {
396  if (procBbOverlaps[proci])
397  {
398  dynSendMap[proci].append(facei);
399  }
400  }
401  }
402  }
403 
404  // Convert dynamicList to labelList
405  sendMap.setSize(Pstream::nProcs());
406  forAll(sendMap, proci)
407  {
408  sendMap[proci].transfer(dynSendMap[proci]);
409  }
410  }
411 
412  // Debug printing
413  if (debug)
414  {
415  Pout<< "Of my " << faces.size() << " I need to send to:" << nl
416  << "\tproc\tfaces" << endl;
417  forAll(sendMap, proci)
418  {
419  Pout<< '\t' << proci << '\t' << sendMap[proci].size() << endl;
420  }
421  }
422 
423 
424  // Send over how many faces I need to receive
425  labelListList sendSizes(Pstream::nProcs());
426  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
427  forAll(sendMap, proci)
428  {
429  sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
430  }
431  Pstream::gatherList(sendSizes);
432  Pstream::scatterList(sendSizes);
433 
434 
435  // Determine order of receiving
436  labelListList constructMap(Pstream::nProcs());
437 
438  // My local segment first
439  constructMap[Pstream::myProcNo()] = identity
440  (
441  sendMap[Pstream::myProcNo()].size()
442  );
443 
444  label segmentI = constructMap[Pstream::myProcNo()].size();
445  forAll(constructMap, proci)
446  {
447  if (proci != Pstream::myProcNo())
448  {
449  // What I need to receive is what other processor is sending to me
450  label nRecv = sendSizes[proci][Pstream::myProcNo()];
451  constructMap[proci].setSize(nRecv);
452 
453  for (label i = 0; i < nRecv; i++)
454  {
455  constructMap[proci][i] = segmentI++;
456  }
457  }
458  }
459 
460  autoPtr<mapDistribute> mapPtr
461  (
462  new mapDistribute
463  (
464  segmentI, // size after construction
465  move(sendMap),
466  move(constructMap)
467  )
468  );
469 
470  return mapPtr;
471 }
472 
473 
474 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
List< treeBoundBox > treeBoundBoxList
List of bounding boxes.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
label nPoints
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const char nl
Definition: Ostream.H:265
Merge points. See below.
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
messageStream Info
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
#define InfoInFunction
Report an information message using Foam::Info.