AMIInterpolationParallelOps.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) 2011-2012 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 template<class SourcePatch, class TargetPatch>
34 (
35  const SourcePatch& srcPatch,
36  const TargetPatch& tgtPatch
37 ) const
38 {
39  label procI = 0;
40 
41  if (Pstream::parRun())
42  {
43  List<label> facesPresentOnProc(Pstream::nProcs(), 0);
44  if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
45  {
46  facesPresentOnProc[Pstream::myProcNo()] = 1;
47  }
48  else
49  {
50  facesPresentOnProc[Pstream::myProcNo()] = 0;
51  }
52 
53  Pstream::gatherList(facesPresentOnProc);
54  Pstream::scatterList(facesPresentOnProc);
55 
56  label nHaveFaces = sum(facesPresentOnProc);
57 
58  if (nHaveFaces > 1)
59  {
60  procI = -1;
61  if (debug)
62  {
63  Info<< "AMIInterpolation::calcDistribution: "
64  << "AMI split across multiple processors" << endl;
65  }
66  }
67  else if (nHaveFaces == 1)
68  {
69  procI = findIndex(facesPresentOnProc, 1);
70  if (debug)
71  {
72  Info<< "AMIInterpolation::calcDistribution: "
73  << "AMI local to processor" << procI << endl;
74  }
75  }
76  }
77 
78 
79  // Either not parallel or no faces on any processor
80  return procI;
81 }
82 
83 
84 template<class SourcePatch, class TargetPatch>
87 (
88  const List<treeBoundBoxList>& procBb,
89  const treeBoundBox& bb,
90  boolList& overlaps
91 ) const
92 {
93  overlaps.setSize(procBb.size());
94  overlaps = false;
95 
96  label nOverlaps = 0;
97 
98  forAll(procBb, procI)
99  {
100  const List<treeBoundBox>& bbs = procBb[procI];
101 
102  forAll(bbs, bbI)
103  {
104  if (bbs[bbI].overlaps(bb))
105  {
106  overlaps[procI] = true;
107  nOverlaps++;
108  break;
109  }
110  }
111  }
112  return nOverlaps;
113 }
114 
115 
116 template<class SourcePatch, class TargetPatch>
118 (
119  const mapDistribute& map,
120  const TargetPatch& pp,
121  const globalIndex& gi,
122  List<faceList>& faces,
123  List<pointField>& points,
124  List<labelList>& faceIDs
125 ) const
126 {
127  PstreamBuffers pBufs(Pstream::nonBlocking);
128 
129  for (label domain = 0; domain < Pstream::nProcs(); domain++)
130  {
131  const labelList& sendElems = map.subMap()[domain];
132 
133  if (domain != Pstream::myProcNo() && sendElems.size())
134  {
135  labelList globalElems(sendElems.size());
136  forAll(sendElems, i)
137  {
138  globalElems[i] = gi.toGlobal(sendElems[i]);
139  }
140 
141  faceList subFaces(UIndirectList<face>(pp, sendElems));
142  primitivePatch subPatch
143  (
144  SubList<face>(subFaces, subFaces.size()),
145  pp.points()
146  );
147 
148  if (debug & 2)
149  {
150  Pout<< "distributePatches: to processor " << domain
151  << " sending faces " << subPatch.faceCentres() << endl;
152  }
153 
154  UOPstream toDomain(domain, pBufs);
155  toDomain
156  << subPatch.localFaces() << subPatch.localPoints()
157  << globalElems;
158  }
159  }
160 
161  // Start receiving
162  pBufs.finishedSends();
163 
164  faces.setSize(Pstream::nProcs());
165  points.setSize(Pstream::nProcs());
166  faceIDs.setSize(Pstream::nProcs());
167 
168  {
169  // Set up 'send' to myself
170  const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
171  faceList subFaces(UIndirectList<face>(pp, sendElems));
172  primitivePatch subPatch
173  (
174  SubList<face>(subFaces, subFaces.size()),
175  pp.points()
176  );
177 
178  // Receive
179  if (debug & 2)
180  {
181  Pout<< "distributePatches: to processor " << Pstream::myProcNo()
182  << " sending faces " << subPatch.faceCentres() << endl;
183  }
184 
185  faces[Pstream::myProcNo()] = subPatch.localFaces();
186  points[Pstream::myProcNo()] = subPatch.localPoints();
187 
188  faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
189  forAll(sendElems, i)
190  {
191  faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
192  }
193  }
194 
195  // Consume
196  for (label domain = 0; domain < Pstream::nProcs(); domain++)
197  {
198  const labelList& recvElems = map.constructMap()[domain];
199 
200  if (domain != Pstream::myProcNo() && recvElems.size())
201  {
202  UIPstream str(domain, pBufs);
203 
204  str >> faces[domain]
205  >> points[domain]
206  >> faceIDs[domain];
207  }
208  }
209 }
210 
211 
212 template<class SourcePatch, class TargetPatch>
215 (
216  const mapDistribute& map,
217  const TargetPatch& tgtPatch,
218  const globalIndex& gi,
219  faceList& tgtFaces,
220  pointField& tgtPoints,
221  labelList& tgtFaceIDs
222 ) const
223 {
224  // Exchange per-processor data
225  List<faceList> allFaces;
226  List<pointField> allPoints;
227  List<labelList> allTgtFaceIDs;
228  distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
229 
230  // Renumber and flatten
231  label nFaces = 0;
232  label nPoints = 0;
233  forAll(allFaces, procI)
234  {
235  nFaces += allFaces[procI].size();
236  nPoints += allPoints[procI].size();
237  }
238 
239  tgtFaces.setSize(nFaces);
240  tgtPoints.setSize(nPoints);
241  tgtFaceIDs.setSize(nFaces);
242 
243  nFaces = 0;
244  nPoints = 0;
245 
246  // My own data first
247  {
248  const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
249  SubList<label>(tgtFaceIDs, faceIDs.size()).assign(faceIDs);
250 
251  const faceList& fcs = allFaces[Pstream::myProcNo()];
252  forAll(fcs, i)
253  {
254  const face& f = fcs[i];
255  face& newF = tgtFaces[nFaces++];
256  newF.setSize(f.size());
257  forAll(f, fp)
258  {
259  newF[fp] = f[fp] + nPoints;
260  }
261  }
262 
263  const pointField& pts = allPoints[Pstream::myProcNo()];
264  forAll(pts, i)
265  {
266  tgtPoints[nPoints++] = pts[i];
267  }
268  }
269 
270 
271  // Other proc data follows
272  forAll(allFaces, procI)
273  {
274  if (procI != Pstream::myProcNo())
275  {
276  const labelList& faceIDs = allTgtFaceIDs[procI];
277  SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces).assign(faceIDs);
278 
279  const faceList& fcs = allFaces[procI];
280  forAll(fcs, i)
281  {
282  const face& f = fcs[i];
283  face& newF = tgtFaces[nFaces++];
284  newF.setSize(f.size());
285  forAll(f, fp)
286  {
287  newF[fp] = f[fp] + nPoints;
288  }
289  }
290 
291  const pointField& pts = allPoints[procI];
292  forAll(pts, i)
293  {
294  tgtPoints[nPoints++] = pts[i];
295  }
296  }
297  }
298 
299  // Merge
300  labelList oldToNew;
301  pointField newTgtPoints;
302  bool hasMerged = mergePoints
303  (
304  tgtPoints,
305  SMALL,
306  false,
307  oldToNew,
308  newTgtPoints
309  );
310 
311  if (hasMerged)
312  {
313  if (debug)
314  {
315  Pout<< "Merged from " << tgtPoints.size()
316  << " down to " << newTgtPoints.size() << " points" << endl;
317  }
318 
319  tgtPoints.transfer(newTgtPoints);
320  forAll(tgtFaces, i)
321  {
322  inplaceRenumber(oldToNew, tgtFaces[i]);
323  }
324  }
325 }
326 
327 
328 template<class SourcePatch, class TargetPatch>
331 (
332  const SourcePatch& srcPatch,
333  const TargetPatch& tgtPatch
334 ) const
335 {
336  // Get decomposition of patch
337  List<treeBoundBoxList> procBb(Pstream::nProcs());
338 
339  if (srcPatch.size())
340  {
341  procBb[Pstream::myProcNo()] = treeBoundBoxList
342  (
343  1, // For now single bounding box per proc
344  treeBoundBox
345  (
346  srcPatch.points(),
347  srcPatch.meshPoints()
348  )
349  );
350  }
351  else
352  {
353  procBb[Pstream::myProcNo()] = treeBoundBoxList();
354  }
355 
356  // slightly increase size of bounding boxes to allow for cases where
357  // bounding boxes are perfectly alligned
358  forAll(procBb[Pstream::myProcNo()], bbI)
359  {
360  treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
361  bb.inflate(0.01);
362  }
363 
364  Pstream::gatherList(procBb);
365  Pstream::scatterList(procBb);
366 
367 
368  if (debug)
369  {
370  Info<< "Determining extent of srcPatch per processor:" << nl
371  << "\tproc\tbb" << endl;
372  forAll(procBb, procI)
373  {
374  Info<< '\t' << procI << '\t' << procBb[procI] << endl;
375  }
376  }
377 
378 
379  // Determine which faces of tgtPatch overlaps srcPatch per proc
380  const faceList& faces = tgtPatch.localFaces();
381  const pointField& points = tgtPatch.localPoints();
382 
383  labelListList sendMap;
384 
385  {
386  // Per processor indices into all segments to send
387  List<DynamicList<label> > dynSendMap(Pstream::nProcs());
388 
389  // Work array - whether processor bb overlaps the face bounds
390  boolList procBbOverlaps(Pstream::nProcs());
391 
392  forAll(faces, faceI)
393  {
394  if (faces[faceI].size())
395  {
396  treeBoundBox faceBb(points, faces[faceI]);
397 
398  // Find the processor this face overlaps
399  calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
400 
401  forAll(procBbOverlaps, procI)
402  {
403  if (procBbOverlaps[procI])
404  {
405  dynSendMap[procI].append(faceI);
406  }
407  }
408  }
409  }
410 
411  // Convert dynamicList to labelList
412  sendMap.setSize(Pstream::nProcs());
413  forAll(sendMap, procI)
414  {
415  sendMap[procI].transfer(dynSendMap[procI]);
416  }
417  }
418 
419  // Debug printing
420  if (debug)
421  {
422  Pout<< "Of my " << faces.size() << " I need to send to:" << nl
423  << "\tproc\tfaces" << endl;
424  forAll(sendMap, procI)
425  {
426  Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
427  }
428  }
429 
430 
431  // Send over how many faces I need to receive
432  labelListList sendSizes(Pstream::nProcs());
433  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
434  forAll(sendMap, procI)
435  {
436  sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
437  }
438  Pstream::gatherList(sendSizes);
439  Pstream::scatterList(sendSizes);
440 
441 
442  // Determine order of receiving
443  labelListList constructMap(Pstream::nProcs());
444 
445  // My local segment first
446  constructMap[Pstream::myProcNo()] = identity
447  (
448  sendMap[Pstream::myProcNo()].size()
449  );
450 
451  label segmentI = constructMap[Pstream::myProcNo()].size();
452  forAll(constructMap, procI)
453  {
454  if (procI != Pstream::myProcNo())
455  {
456  // What I need to receive is what other processor is sending to me
457  label nRecv = sendSizes[procI][Pstream::myProcNo()];
458  constructMap[procI].setSize(nRecv);
459 
460  for (label i = 0; i < nRecv; i++)
461  {
462  constructMap[procI][i] = segmentI++;
463  }
464  }
465  }
466 
467  autoPtr<mapDistribute> mapPtr
468  (
469  new mapDistribute
470  (
471  segmentI, // size after construction
472  sendMap.xfer(),
473  constructMap.xfer()
474  )
475  );
476 
477  return mapPtr;
478 }
479 
480 
481 // ************************************************************************* //
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
labelList f(nPoints)
Merge points. See below.
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
List< face > faceList
Definition: faceListFwd.H:43
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
List< treeBoundBox > treeBoundBoxList
List of bounding boxes.
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
#define forAll(list, i)
Definition: UList.H:421
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)
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
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.
List< label > labelList
A List of labels.
Definition: labelList.H:56
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
label nPoints
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50