FreeFOAM The Cross-Platform CFD Toolkit
globalPoints.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) 1991-2010 OpenCFD Ltd.
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 "globalPoints.H"
29 #include <OpenFOAM/polyMesh.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
34 
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 // Total number of points on processor patches. Is upper limit for number
39 // of shared points
40 Foam::label Foam::globalPoints::countPatchPoints
41 (
42  const polyBoundaryMesh& patches
43 )
44 {
45  label nTotPoints = 0;
46 
47  forAll(patches, patchI)
48  {
49  const polyPatch& pp = patches[patchI];
50 
51  if
52  (
53  (Pstream::parRun() && isA<processorPolyPatch>(pp))
54  || isA<cyclicPolyPatch>(pp)
55  )
56  {
57  nTotPoints += pp.nPoints();
58  }
59  }
60 
61  return nTotPoints;
62 }
63 
64 
65 // Collect all topological information about a point on a patch.
66 // (this information is the patch faces using the point and the relative
67 // position of the point in the face)
68 void Foam::globalPoints::addToSend
69 (
70  const primitivePatch& pp,
71  const label patchPointI,
72  const procPointList& knownInfo,
73 
74  DynamicList<label>& patchFaces,
75  DynamicList<label>& indexInFace,
76  DynamicList<procPointList>& allInfo
77 )
78 {
79  label meshPointI = pp.meshPoints()[patchPointI];
80 
81  // Add all faces using the point so we are sure we find it on the
82  // other side.
83  const labelList& pFaces = pp.pointFaces()[patchPointI];
84 
85  forAll(pFaces, i)
86  {
87  label patchFaceI = pFaces[i];
88 
89  const face& f = pp[patchFaceI];
90 
91  patchFaces.append(patchFaceI);
92  indexInFace.append(findIndex(f, meshPointI));
93  allInfo.append(knownInfo);
94  }
95 }
96 
97 
98 // Add nbrInfo to myInfo. Return true if anything changed.
99 // nbrInfo is for a point a list of all the processors using it (and the
100 // meshPoint label on that processor)
101 bool Foam::globalPoints::mergeInfo
102 (
103  const procPointList& nbrInfo,
104  procPointList& myInfo
105 )
106 {
107  // Indices of entries in nbrInfo not yet in myInfo.
108  DynamicList<label> newInfo(nbrInfo.size());
109 
110  forAll(nbrInfo, i)
111  {
112  const procPoint& info = nbrInfo[i];
113 
114  // Check if info already in myInfo.
115  label index = -1;
116 
117  forAll(myInfo, j)
118  {
119  if (myInfo[j] == info)
120  {
121  // Already have information for processor/point combination
122  // in my list so skip.
123  index = j;
124 
125  break;
126  }
127  }
128 
129  if (index == -1)
130  {
131  // Mark this information as being not yet in myInfo
132  newInfo.append(i);
133  }
134  }
135 
136  newInfo.shrink();
137 
138  // Append all nbrInfos referenced in newInfo to myInfo.
139 
140  label index = myInfo.size();
141 
142  myInfo.setSize(index + newInfo.size());
143 
144  forAll(newInfo, i)
145  {
146  myInfo[index++] = nbrInfo[newInfo[i]];
147  }
148 
149  // Did anything change?
150  return newInfo.size() > 0;
151 }
152 
153 
154 // Updates database of current information on meshpoints with nbrInfo.
155 // Uses mergeInfo above. Returns true if data kept for meshPointI changed.
156 bool Foam::globalPoints::storeInfo
157 (
158  const procPointList& nbrInfo,
159  const label meshPointI
160 )
161 {
162  label infoChanged = false;
163 
164  // Get the index into the procPoints list.
165  Map<label>::iterator iter = meshToProcPoint_.find(meshPointI);
166 
167  if (iter != meshToProcPoint_.end())
168  {
169  procPointList& knownInfo = procPoints_[iter()];
170 
171  if (mergeInfo(nbrInfo, knownInfo))
172  {
173  infoChanged = true;
174  }
175  }
176  else
177  {
178  procPointList knownInfo(1);
179  knownInfo[0][0] = Pstream::myProcNo();
180  knownInfo[0][1] = meshPointI;
181 
182  if (mergeInfo(nbrInfo, knownInfo))
183  {
184  // Update addressing from into procPoints
185  meshToProcPoint_.insert(meshPointI, procPoints_.size());
186  // Insert into list of equivalences.
187  procPoints_.append(knownInfo);
188 
189  infoChanged = true;
190  }
191  }
192  return infoChanged;
193 }
194 
195 
196 // Insert my own points into structure and mark as changed.
197 void Foam::globalPoints::initOwnPoints
198 (
199  const bool allPoints,
200  labelHashSet& changedPoints
201 )
202 {
203  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
204 
205  forAll(patches, patchI)
206  {
207  const polyPatch& pp = patches[patchI];
208 
209  if
210  (
211  (Pstream::parRun() && isA<processorPolyPatch>(pp))
212  || isA<cyclicPolyPatch>(pp)
213  )
214  {
215  const labelList& meshPoints = pp.meshPoints();
216 
217  if (allPoints)
218  {
219  // All points on patch
220  forAll(meshPoints, i)
221  {
222  label meshPointI = meshPoints[i];
223 
224  procPointList knownInfo(1);
225  knownInfo[0][0] = Pstream::myProcNo();
226  knownInfo[0][1] = meshPointI;
227 
228  // Update addressing from meshpoint to index in procPoints
229  meshToProcPoint_.insert(meshPointI, procPoints_.size());
230  // Insert into list of equivalences.
231  procPoints_.append(knownInfo);
232 
233  // Update changedpoints info.
234  changedPoints.insert(meshPointI);
235  }
236  }
237  else
238  {
239  // Boundary points only
240  const labelList& boundaryPoints = pp.boundaryPoints();
241 
242  forAll(boundaryPoints, i)
243  {
244  label meshPointI = meshPoints[boundaryPoints[i]];
245 
246  procPointList knownInfo(1);
247  knownInfo[0][0] = Pstream::myProcNo();
248  knownInfo[0][1] = meshPointI;
249 
250  // Update addressing from meshpoint to index in procPoints
251  meshToProcPoint_.insert(meshPointI, procPoints_.size());
252  // Insert into list of equivalences.
253  procPoints_.append(knownInfo);
254 
255  // Update changedpoints info.
256  changedPoints.insert(meshPointI);
257  }
258  }
259  }
260  }
261 }
262 
263 
264 // Send all my info on changedPoints_ to my neighbours.
265 void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
266  const
267 {
268  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
269 
270  forAll(patches, patchI)
271  {
272  const polyPatch& pp = patches[patchI];
273 
274  if (Pstream::parRun() && isA<processorPolyPatch>(pp))
275  {
276  // Information to send:
277  // patch face
278  DynamicList<label> patchFaces(pp.nPoints());
279  // index in patch face
280  DynamicList<label> indexInFace(pp.nPoints());
281  // all information I currently hold about this patchPoint
282  DynamicList<procPointList> allInfo(pp.nPoints());
283 
284 
285  // Now collect information on all mesh points mentioned in
286  // changedPoints. Note that these points only should occur on
287  // processorPatches (or rather this is a limitation!).
288 
289  const labelList& meshPoints = pp.meshPoints();
290 
291  forAll(meshPoints, patchPointI)
292  {
293  label meshPointI = meshPoints[patchPointI];
294 
295  if (changedPoints.found(meshPointI))
296  {
297  label index = meshToProcPoint_[meshPointI];
298 
299  const procPointList& knownInfo = procPoints_[index];
300 
301  // Add my information about meshPointI to the send buffers
302  addToSend
303  (
304  pp,
305  patchPointI,
306  knownInfo,
307 
308  patchFaces,
309  indexInFace,
310  allInfo
311  );
312  }
313  }
314  patchFaces.shrink();
315  indexInFace.shrink();
316  allInfo.shrink();
317 
318  // Send to neighbour
319  {
320  const processorPolyPatch& procPatch =
321  refCast<const processorPolyPatch>(pp);
322 
323  if (debug)
324  {
325  Pout<< " Sending to "
326  << procPatch.neighbProcNo() << " point information:"
327  << patchFaces.size() << endl;
328  }
329 
330  OPstream toNeighbour
331  (
333  procPatch.neighbProcNo()
334  );
335 
336  toNeighbour << patchFaces << indexInFace << allInfo;
337  }
338  }
339  }
340 }
341 
342 
343 // Receive all my neighbours' information and merge with mine.
344 // After finishing will have updated
345 // - procPoints_ : all neighbour information merged in.
346 // - meshToProcPoint_
347 // - changedPoints: all meshPoints for which something changed.
348 void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
349 {
350  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
351 
352  // Reset changed points
353  changedPoints.clear();
354 
355  forAll(patches, patchI)
356  {
357  const polyPatch& pp = patches[patchI];
358 
359  if (Pstream::parRun() && isA<processorPolyPatch>(pp))
360  {
361  const processorPolyPatch& procPatch =
362  refCast<const processorPolyPatch>(pp);
363 
364  labelList patchFaces;
365  labelList indexInFace;
366  List<procPointList> nbrInfo;
367 
368  {
369  IPstream fromNeighbour
370  (
372  procPatch.neighbProcNo()
373  );
374  fromNeighbour >> patchFaces >> indexInFace >> nbrInfo;
375  }
376 
377  if (debug)
378  {
379  Pout<< " Received from "
380  << procPatch.neighbProcNo() << " point information:"
381  << patchFaces.size() << endl;
382  }
383 
384  forAll(patchFaces, i)
385  {
386  const face& f = pp[patchFaces[i]];
387 
388  // Get index in this face from index on face on other side.
389  label index = (f.size() - indexInFace[i]) % f.size();
390 
391  // Get the meshpoint on my side
392  label meshPointI = f[index];
393 
394  //const procPointList& info = nbrInfo[i];
395  //Pout << "Received for my coord "
396  // << mesh_.points()[meshPointI];
397  //
398  //forAll(info, j)
399  //{
400  // Pout<< ' ' <<info[j];
401  //}
402  //Pout<< endl;
403 
404  if (storeInfo(nbrInfo[i], meshPointI))
405  {
406  changedPoints.insert(meshPointI);
407  }
408  }
409  }
410  else if (isA<cyclicPolyPatch>(pp))
411  {
412  // Handle cyclics: send lower half to upper half and vice versa.
413  // Or since they both are in memory just do it point by point.
414 
415  const cyclicPolyPatch& cycPatch =
416  refCast<const cyclicPolyPatch>(pp);
417 
418  const labelList& meshPoints = pp.meshPoints();
419 
420  //const edgeList& connections = cycPatch.coupledPoints();
421  const edgeList connections(coupledPoints(cycPatch));
422 
423  forAll(connections, i)
424  {
425  const edge& e = connections[i];
426 
427  label meshPointA = meshPoints[e[0]];
428  label meshPointB = meshPoints[e[1]];
429 
430  // Do we have information on meshPointA?
431  Map<label>::iterator procPointA =
432  meshToProcPoint_.find(meshPointA);
433 
434  if (procPointA != meshToProcPoint_.end())
435  {
436  // Store A info onto pointB
437  if (storeInfo(procPoints_[procPointA()], meshPointB))
438  {
439  changedPoints.insert(meshPointB);
440  }
441  }
442 
443  // Same for info on pointB
444  Map<label>::iterator procPointB =
445  meshToProcPoint_.find(meshPointB);
446 
447  if (procPointB != meshToProcPoint_.end())
448  {
449  // Store B info onto pointA
450  if (storeInfo(procPoints_[procPointB()], meshPointA))
451  {
452  changedPoints.insert(meshPointA);
453  }
454  }
455  }
456  }
457  }
458 }
459 
460 
461 // Remove entries which are handled by normal face-face communication. I.e.
462 // those points where the equivalence list is only me and my (face)neighbour
463 void Foam::globalPoints::remove(const Map<label>& directNeighbours)
464 {
465  // Save old ones.
466  Map<label> oldMeshToProcPoint(meshToProcPoint_);
467  meshToProcPoint_.clear();
468 
469  List<procPointList> oldProcPoints;
470  oldProcPoints.transfer(procPoints_);
471 
472  // Go through all equivalences
473  for
474  (
475  Map<label>::const_iterator iter = oldMeshToProcPoint.begin();
476  iter != oldMeshToProcPoint.end();
477  ++iter
478  )
479  {
480  label meshPointI = iter.key();
481  const procPointList& pointInfo = oldProcPoints[iter()];
482 
483  if (pointInfo.size() == 2)
484  {
485  // I will be in this equivalence list.
486  // Check whether my direct (=face) neighbour
487  // is in it. This would be an ordinary connection and can be
488  // handled by normal face-face connectivity.
489 
490  const procPoint& a = pointInfo[0];
491  const procPoint& b = pointInfo[1];
492 
493  if
494  (
495  (a[0] == Pstream::myProcNo() && directNeighbours.found(a[1]))
496  || (b[0] == Pstream::myProcNo() && directNeighbours.found(b[1]))
497  )
498  {
499  // Normal faceNeighbours
500  if (a[0] == Pstream::myProcNo())
501  {
502  //Pout<< "Removing direct neighbour:"
503  // << mesh_.points()[a[1]]
504  // << endl;
505  }
506  else if (b[0] == Pstream::myProcNo())
507  {
508  //Pout<< "Removing direct neighbour:"
509  // << mesh_.points()[b[1]]
510  // << endl;
511  }
512  }
513  else
514  {
515  // This condition will be very rare: points are used by
516  // two processors which are not face-face connected.
517  // e.g.
518  // +------+------+
519  // | wall | B |
520  // +------+------+
521  // | A | wall |
522  // +------+------+
523  // Processor A and B share a point. Note that this only will
524  // be found if the two domains are face connected at all
525  // (not shown in the picture)
526 
527  meshToProcPoint_.insert(meshPointI, procPoints_.size());
528  procPoints_.append(pointInfo);
529  }
530  }
531  else if (pointInfo.size() == 1)
532  {
533  // This happens for 'wedge' like cyclics where the two halves
534  // come together in the same point so share the same meshPoint.
535  // So this meshPoint will have info of size one only.
536  if
537  (
538  pointInfo[0][0] != Pstream::myProcNo()
539  || !directNeighbours.found(pointInfo[0][1])
540  )
541  {
542  meshToProcPoint_.insert(meshPointI, procPoints_.size());
543  procPoints_.append(pointInfo);
544  }
545  }
546  else
547  {
548  meshToProcPoint_.insert(meshPointI, procPoints_.size());
549  procPoints_.append(pointInfo);
550  }
551  }
552 
553  procPoints_.shrink();
554 }
555 
556 
557 // Get (indices of) points for which I am master (= lowest numbered point on
558 // lowest numbered processor).
559 // (equivalence lists should be complete by now)
560 Foam::labelList Foam::globalPoints::getMasterPoints() const
561 {
562  labelList masterPoints(nPatchPoints_);
563  label nMaster = 0;
564 
565  // Go through all equivalences and determine meshPoints where I am master.
566  for
567  (
568  Map<label>::const_iterator iter = meshToProcPoint_.begin();
569  iter != meshToProcPoint_.end();
570  ++iter
571  )
572  {
573  label meshPointI = iter.key();
574  const procPointList& pointInfo = procPoints_[iter()];
575 
576  if (pointInfo.size() < 2)
577  {
578  // Points should have an equivalence list >= 2 since otherwise
579  // they would be direct neighbours and have been removed in the
580  // call to 'remove'.
581  Pout<< "MeshPoint:" << meshPointI
582  << " coord:" << mesh_.points()[meshPointI]
583  << " has no corresponding point on a neighbouring processor"
584  << endl;
585  FatalErrorIn("globalPoints::getMasterPoints()")
586  << '[' << Pstream::myProcNo() << ']'
587  << " MeshPoint:" << meshPointI
588  << " coord:" << mesh_.points()[meshPointI]
589  << " has no corresponding point on a neighbouring processor"
590  << abort(FatalError);
591  }
592  else
593  {
594  // Find lowest processor and lowest mesh point on this processor.
595  label lowestProcI = labelMax;
596  label lowestPointI = labelMax;
597 
598  forAll(pointInfo, i)
599  {
600  label proc = pointInfo[i][0];
601 
602  if (proc < lowestProcI)
603  {
604  lowestProcI = proc;
605  lowestPointI = pointInfo[i][1];
606  }
607  else if (proc == lowestProcI)
608  {
609  lowestPointI = min(lowestPointI, pointInfo[i][1]);
610  }
611  }
612 
613  if
614  (
615  lowestProcI == Pstream::myProcNo()
616  && lowestPointI == meshPointI
617  )
618  {
619  // I am lowest numbered processor and point. Add to my list.
620  masterPoints[nMaster++] = meshPointI;
621  }
622  }
623  }
624 
625  masterPoints.setSize(nMaster);
626 
627  return masterPoints;
628 }
629 
630 
631 // Send subset of lists
632 void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
633 {
634  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
635 
636  forAll(patches, patchI)
637  {
638  const polyPatch& pp = patches[patchI];
639 
640  if (Pstream::parRun() && isA<processorPolyPatch>(pp))
641  {
642  const processorPolyPatch& procPatch =
643  refCast<const processorPolyPatch>(pp);
644 
645  OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
646 
647  if (debug)
648  {
649  Pout<< "Sending to " << procPatch.neighbProcNo()
650  << " changed sharedPoints info:"
651  << changedIndices.size() << endl;
652  }
653 
654  toNeighbour
655  << UIndirectList<label>(sharedPointAddr_, changedIndices)()
656  << UIndirectList<label>(sharedPointLabels_, changedIndices)();
657  }
658  }
659 }
660 
661 
662 // Receive shared point indices for all my shared points. Note that since
663 // there are only a few here we can build a reverse map using the meshpoint
664 // instead of doing all this relative point indexing (patch face + index in
665 // face) as in send/receivePatchPoints
666 void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
667 {
668  changedIndices.setSize(sharedPointAddr_.size());
669  label nChanged = 0;
670 
671  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
672 
673  // Receive and set shared points
674  forAll(patches, patchI)
675  {
676  const polyPatch& pp = patches[patchI];
677 
678  if (Pstream::parRun() && isA<processorPolyPatch>(pp))
679  {
680  const processorPolyPatch& procPatch =
681  refCast<const processorPolyPatch>(pp);
682 
683  // Map from neighbouring meshPoint to sharedPoint)
684  Map<label> nbrSharedPoints(sharedPointAddr_.size());
685 
686  {
687  // Receive meshPoints on neighbour and sharedPoints and build
688  // map from it. Note that we could have built the map on the
689  // neighbour and sent it over.
690  labelList nbrSharedPointAddr;
691  labelList nbrSharedPointLabels;
692 
693  {
694  IPstream fromNeighbour
695  (
697  procPatch.neighbProcNo()
698  );
699  fromNeighbour >> nbrSharedPointAddr >> nbrSharedPointLabels;
700  }
701 
702  // Insert into to map
703  forAll(nbrSharedPointLabels, i)
704  {
705  nbrSharedPoints.insert
706  (
707  nbrSharedPointLabels[i], // meshpoint on neighbour
708  nbrSharedPointAddr[i] // sharedPoint label
709  );
710  }
711  }
712 
713 
714  // Merge into whatever information I hold.
715  for
716  (
717  Map<label>::const_iterator iter = meshToProcPoint_.begin();
718  iter != meshToProcPoint_.end();
719  ++iter
720  )
721  {
722  label meshPointI = iter.key();
723  label index = iter();
724 
725  if (sharedPointAddr_[index] == -1)
726  {
727  // No shared point known yet for this meshPoint.
728  // See if was received from neighbour.
729  const procPointList& knownInfo = procPoints_[index];
730 
731  // Check through the whole equivalence list for any
732  // point from the neighbour.
733  forAll(knownInfo, j)
734  {
735  const procPoint& info = knownInfo[j];
736 
737  if
738  (
739  (info[0] == procPatch.neighbProcNo())
740  && nbrSharedPoints.found(info[1])
741  )
742  {
743  // So this knownInfo contains the neighbour point
744  label sharedPointI = nbrSharedPoints[info[1]];
745 
746  sharedPointAddr_[index] = sharedPointI;
747  sharedPointLabels_[index] = meshPointI;
748  changedIndices[nChanged++] = index;
749 
750  break;
751  }
752  }
753  }
754  }
755  }
756  else if (isA<cyclicPolyPatch>(pp))
757  {
758  const cyclicPolyPatch& cycPatch =
759  refCast<const cyclicPolyPatch>(pp);
760 
761  // Build map from meshPoint to sharedPoint
762  Map<label> meshToSharedPoint(sharedPointAddr_.size());
763  forAll(sharedPointLabels_, i)
764  {
765  label meshPointI = sharedPointLabels_[i];
766 
767  meshToSharedPoint.insert(meshPointI, sharedPointAddr_[i]);
768  }
769 
770  // Sync all info.
771  //const edgeList& connections = cycPatch.coupledPoints();
772  const edgeList connections(coupledPoints(cycPatch));
773 
774  forAll(connections, i)
775  {
776  const edge& e = connections[i];
777 
778  label meshPointA = pp.meshPoints()[e[0]];
779  label meshPointB = pp.meshPoints()[e[1]];
780 
781  // Do we already have shared point for meshPointA?
782  Map<label>::iterator fndA = meshToSharedPoint.find(meshPointA);
783  Map<label>::iterator fndB = meshToSharedPoint.find(meshPointB);
784 
785  if (fndA != meshToSharedPoint.end())
786  {
787  if (fndB != meshToSharedPoint.end())
788  {
789  if (fndA() != fndB())
790  {
792  (
793  "globalPoints::receiveSharedPoints"
794  "(labelList&)"
795  ) << "On patch " << pp.name()
796  << " connected points " << meshPointA
797  << ' ' << mesh_.points()[meshPointA]
798  << " and " << meshPointB
799  << ' ' << mesh_.points()[meshPointB]
800  << " are mapped to different shared points: "
801  << fndA() << " and " << fndB()
802  << abort(FatalError);
803  }
804  }
805  else
806  {
807  // No shared point yet for B.
808  label sharedPointI = fndA();
809 
810  // Store shared point for meshPointB
811  label index = meshToProcPoint_[meshPointB];
812 
813  sharedPointAddr_[index] = sharedPointI;
814  sharedPointLabels_[index] = meshPointB;
815  changedIndices[nChanged++] = index;
816  }
817  }
818  else
819  {
820  // No shared point yet for A.
821  if (fndB != meshToSharedPoint.end())
822  {
823  label sharedPointI = fndB();
824 
825  // Store shared point for meshPointA
826  label index = meshToProcPoint_[meshPointA];
827 
828  sharedPointAddr_[index] = sharedPointI;
829  sharedPointLabels_[index] = meshPointA;
830  changedIndices[nChanged++] = index;
831  }
832  }
833  }
834  }
835  }
836 
837  changedIndices.setSize(nChanged);
838 }
839 
840 
841 Foam::edgeList Foam::globalPoints::coupledPoints(const cyclicPolyPatch& pp)
842 {
843  // Look at cyclic patch as two halves, A and B.
844  // Now all we know is that relative face index in halfA is same
845  // as coupled face in halfB and also that the 0th vertex
846  // corresponds.
847 
848  // From halfA point to halfB or -1.
849  labelList coupledPoint(pp.nPoints(), -1);
850 
851  for (label patchFaceA = 0; patchFaceA < pp.size()/2; patchFaceA++)
852  {
853  const face& fA = pp.localFaces()[patchFaceA];
854 
855  forAll(fA, indexA)
856  {
857  label patchPointA = fA[indexA];
858 
859  if (coupledPoint[patchPointA] == -1)
860  {
861  const face& fB = pp.localFaces()[patchFaceA + pp.size()/2];
862 
863  label indexB = (fB.size() - indexA) % fB.size();
864 
865  coupledPoint[patchPointA] = fB[indexB];
866  }
867  }
868  }
869 
870  edgeList connected(pp.nPoints());
871 
872  // Extract coupled points.
873  label connectedI = 0;
874 
875  forAll(coupledPoint, i)
876  {
877  if (coupledPoint[i] != -1)
878  {
879  connected[connectedI++] = edge(i, coupledPoint[i]);
880  }
881  }
882 
883  connected.setSize(connectedI);
884 
885  return connected;
886 }
887 
888 
889 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
890 
891 // Construct from components
892 Foam::globalPoints::globalPoints(const polyMesh& mesh)
893 :
894  mesh_(mesh),
895  nPatchPoints_(countPatchPoints(mesh.boundaryMesh())),
896  procPoints_(nPatchPoints_),
897  meshToProcPoint_(nPatchPoints_),
898  sharedPointAddr_(0),
899  sharedPointLabels_(0),
900  nGlobalPoints_(0)
901 {
902  if (debug)
903  {
904  Pout<< "globalPoints::globalPoints(const polyMesh&) : "
905  << "doing processor to processor communication to get sharedPoints"
906  << endl;
907  }
908 
909  labelHashSet changedPoints(nPatchPoints_);
910 
911  // Initialize procPoints with my patch points. Keep track of points
912  // inserted (in changedPoints)
913  // There are two possible forms of this:
914  // - initialize with all patch points (allPoints = true). This causes all
915  // patch points to be exchanged so a lot of information gets stored and
916  // transferred. This all gets filtered out later when removing the
917  // equivalence lists of size 2.
918  // - initialize with boundary points of patches only (allPoints = false).
919  // This should work for all decompositions except extreme ones where a
920  // shared point is not on the boundary of any processor patches using it.
921  // This would happen if a domain was pinched such that two patches share
922  // a point or edge.
923  initOwnPoints(true, changedPoints);
924 
925  // Do one exchange iteration to get neighbour points.
926  sendPatchPoints(changedPoints);
927  receivePatchPoints(changedPoints);
928 
929 
930  // Save neighbours reachable through face-face communication.
931  Map<label> neighbourList(meshToProcPoint_);
932 
933 
934  // Exchange until nothing changes on all processors.
935  bool changed = false;
936 
937  do
938  {
939  sendPatchPoints(changedPoints);
940  receivePatchPoints(changedPoints);
941 
942  changed = changedPoints.size() > 0;
943  reduce(changed, orOp<bool>());
944 
945  } while (changed);
946 
947 
948  // Remove direct neighbours from point equivalences.
949  remove(neighbourList);
950 
951 
952  //Pout<< "After removing locally connected points:" << endl;
953  //for
954  //(
955  // Map<label>::const_iterator iter = meshToProcPoint_.begin();
956  // iter != meshToProcPoint_.end();
957  // ++iter
958  //)
959  //{
960  // label meshPointI = iter.key();
961  // const procPointList& pointInfo = procPoints_[iter()];
962  //
963  // forAll(pointInfo, i)
964  // {
965  // Pout<< " pointI:" << meshPointI << ' '
966  // << mesh.points()[meshPointI]
967  // << " connected to proc " << pointInfo[i][0]
968  // << " point:" << pointInfo[i][1]
969  // << endl;
970  // }
971  //}
972 
973 
974  // We now have - in procPoints_ - a list of points which are shared between
975  // multiple processors. These are the ones for which are sharedPoint
976  // needs to be determined. This is done by having the lowest numbered
977  // processor in the equivalence list 'ask' for a sharedPoint number
978  // and then distribute it across processor patches to the non-master
979  // processors. Note: below piece of coding is not very efficient. Uses
980  // a Map where possibly it shouldn't
981 
982  // Initialize sharedPoint addressing. Is for every entry in procPoints_
983  // the sharedPoint.
984  sharedPointAddr_.setSize(meshToProcPoint_.size());
985  sharedPointAddr_ = -1;
986  sharedPointLabels_.setSize(meshToProcPoint_.size());
987  sharedPointLabels_ = -1;
988 
989 
990  // Get points for which I am master (lowest numbered proc)
991  labelList masterPoints(getMasterPoints());
992 
993 
994  // Determine number of master points on all processors.
995  labelList sharedPointSizes(Pstream::nProcs());
996  sharedPointSizes[Pstream::myProcNo()] = masterPoints.size();
997 
998  Pstream::gatherList(sharedPointSizes);
999  Pstream::scatterList(sharedPointSizes);
1000 
1001  if (debug)
1002  {
1003  Pout<< "sharedPointSizes:" << sharedPointSizes << endl;
1004  }
1005 
1006  // Get total number of shared points
1007  nGlobalPoints_ = 0;
1008  forAll(sharedPointSizes, procI)
1009  {
1010  nGlobalPoints_ += sharedPointSizes[procI];
1011  }
1012 
1013  // Assign sharedPoint labels. Every processor gets assigned consecutive
1014  // numbers for its master points.
1015  // These are assigned in processor order so processor0 gets
1016  // 0..masterPoints.size()-1 etc.
1017 
1018  // My point labels start after those of lower numbered processors
1019  label sharedPointI = 0;
1020  for (label procI = 0; procI < Pstream::myProcNo(); procI++)
1021  {
1022  sharedPointI += sharedPointSizes[procI];
1023  }
1024 
1025  forAll(masterPoints, i)
1026  {
1027  label meshPointI = masterPoints[i];
1028 
1029  label index = meshToProcPoint_[meshPointI];
1030 
1031  sharedPointLabels_[index] = meshPointI;
1032  sharedPointAddr_[index] = sharedPointI++;
1033  }
1034 
1035 
1036  // Now we have a sharedPointLabel for some of the entries in procPoints.
1037  // Send this information to neighbours. Receive their information.
1038  // Loop until nothing changes.
1039 
1040  // Initial subset to send is points for which I have sharedPoints
1041  labelList changedIndices(sharedPointAddr_.size());
1042  label nChanged = 0;
1043 
1044  forAll(sharedPointAddr_, i)
1045  {
1046  if (sharedPointAddr_[i] != -1)
1047  {
1048  changedIndices[nChanged++] = i;
1049  }
1050  }
1051  changedIndices.setSize(nChanged);
1052 
1053  changed = false;
1054 
1055  do
1056  {
1057  if (debug)
1058  {
1059  Pout<< "Determined " << changedIndices.size() << " shared points."
1060  << " Exchanging them" << endl;
1061  }
1062  sendSharedPoints(changedIndices);
1063  receiveSharedPoints(changedIndices);
1064 
1065  changed = changedIndices.size() > 0;
1066  reduce(changed, orOp<bool>());
1067 
1068  } while (changed);
1069 
1070 
1071  forAll(sharedPointLabels_, i)
1072  {
1073  if (sharedPointLabels_[i] == -1)
1074  {
1075  FatalErrorIn("globalPoints::globalPoints(const polyMesh& mesh)")
1076  << "Problem: shared point on processor " << Pstream::myProcNo()
1077  << " not set at index " << sharedPointLabels_[i] << endl
1078  << "This might mean the individual processor domains are not"
1079  << " connected and the overall domain consists of multiple"
1080  << " regions. You can check this with checkMesh"
1081  << abort(FatalError);
1082  }
1083  }
1084 
1085  if (debug)
1086  {
1087  Pout<< "globalPoints::globalPoints(const polyMesh&) : "
1088  << "Finished global points" << endl;
1089  }
1090 }
1091 
1092 
1093 // ************************ vim: set sw=4 sts=4 et: ************************ //