globalMeshData.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 "globalMeshData.H"
27 #include "Pstream.H"
29 #include "processorPolyPatch.H"
30 #include "globalPoints.H"
31 #include "polyMesh.H"
32 #include "mapDistribute.H"
33 #include "labelIOList.H"
34 #include "mergePoints.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 defineTypeNameAndDebug(globalMeshData, 0);
42 
43 const scalar globalMeshData::matchTol_ = 1e-8;
44 
45 template<>
46 class minEqOp<labelPair>
47 {
48 public:
49  void operator()(labelPair& x, const labelPair& y) const
50  {
51  x[0] = min(x[0], y[0]);
52  x[1] = min(x[1], y[1]);
53  }
54 };
55 }
56 
57 
58 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
59 
60 void Foam::globalMeshData::initProcAddr()
61 {
62  processorPatchIndices_.setSize(mesh_.boundaryMesh().size());
63  processorPatchIndices_ = -1;
64 
65  processorPatchNeighbours_.setSize(mesh_.boundaryMesh().size());
66  processorPatchNeighbours_ = -1;
67 
68  // Construct processor patch indexing. processorPatchNeighbours_ only
69  // set if running in parallel!
70  processorPatches_.setSize(mesh_.boundaryMesh().size());
71 
72  label nNeighbours = 0;
73 
74  forAll(mesh_.boundaryMesh(), patchi)
75  {
76  if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
77  {
78  processorPatches_[nNeighbours] = patchi;
79  processorPatchIndices_[patchi] = nNeighbours++;
80  }
81  }
82  processorPatches_.setSize(nNeighbours);
83 
84 
85  if (Pstream::parRun())
86  {
88 
89  // Send indices of my processor patches to my neighbours
90  forAll(processorPatches_, i)
91  {
92  label patchi = processorPatches_[i];
93 
94  UOPstream toNeighbour
95  (
96  refCast<const processorPolyPatch>
97  (
98  mesh_.boundaryMesh()[patchi]
99  ).neighbProcNo(),
100  pBufs
101  );
102 
103  toNeighbour << processorPatchIndices_[patchi];
104  }
105 
106  pBufs.finishedSends();
107 
108  forAll(processorPatches_, i)
109  {
110  label patchi = processorPatches_[i];
111 
112  UIPstream fromNeighbour
113  (
114  refCast<const processorPolyPatch>
115  (
116  mesh_.boundaryMesh()[patchi]
117  ).neighbProcNo(),
118  pBufs
119  );
120 
121  fromNeighbour >> processorPatchNeighbours_[patchi];
122  }
123  }
124 }
125 
126 
127 void Foam::globalMeshData::calcSharedPoints() const
128 {
129  if
130  (
131  nGlobalPoints_ != -1
132  || sharedPointLabelsPtr_.valid()
133  || sharedPointAddrPtr_.valid()
134  )
135  {
137  << "Shared point addressing already done" << abort(FatalError);
138  }
139 
140  // Calculate all shared points (exclude points that are only
141  // on two coupled patches). This does all the hard work.
142  globalPoints parallelPoints(mesh_, false, true);
143 
144  // Count the number of master points
145  label nMaster = 0;
146  forAll(parallelPoints.pointPoints(), i)
147  {
148  const labelList& pPoints = parallelPoints.pointPoints()[i];
149  const labelList& transPPoints =
150  parallelPoints.transformedPointPoints()[i];
151 
152  if (pPoints.size()+transPPoints.size() > 0)
153  {
154  nMaster++;
155  }
156  }
157 
158  // Allocate global numbers
159  globalIndex masterNumbering(nMaster);
160 
161  nGlobalPoints_ = masterNumbering.size();
162 
163 
164  // Push master number to slaves
165  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166  // 1. Fill master and slave slots
167  nMaster = 0;
168  labelList master(parallelPoints.map().constructSize(), -1);
169  forAll(parallelPoints.pointPoints(), i)
170  {
171  const labelList& pPoints = parallelPoints.pointPoints()[i];
172  const labelList& transPPoints =
173  parallelPoints.transformedPointPoints()[i];
174 
175  if (pPoints.size()+transPPoints.size() > 0)
176  {
177  master[i] = masterNumbering.toGlobal(nMaster);
178  forAll(pPoints, j)
179  {
180  master[pPoints[j]] = master[i];
181  }
182  forAll(transPPoints, j)
183  {
184  master[transPPoints[j]] = master[i];
185  }
186  nMaster++;
187  }
188  }
189 
190 
191  // 2. Push slave slots back to local storage on originating processor
192  // For all the four types of points:
193  // - local master : already set
194  // - local transformed slave point : the reverse transform at
195  // reverseDistribute will have copied it back to its originating local
196  // point
197  // - remote untransformed slave point : sent back to originating processor
198  // - remote transformed slave point : the reverse transform will
199  // copy it back into the remote slot which then gets sent back to
200  // originating processor
201 
202  parallelPoints.map().reverseDistribute
203  (
204  parallelPoints.map().constructSize(),
205  master
206  );
207 
208 
209  // Collect all points that are a master or refer to a master.
210  nMaster = 0;
211  forAll(parallelPoints.pointPoints(), i)
212  {
213  if (master[i] != -1)
214  {
215  nMaster++;
216  }
217  }
218 
219  sharedPointLabelsPtr_.reset(new labelList(nMaster));
220  labelList& sharedPointLabels = sharedPointLabelsPtr_();
221  sharedPointAddrPtr_.reset(new labelList(nMaster));
222  labelList& sharedPointAddr = sharedPointAddrPtr_();
223  nMaster = 0;
224 
225  forAll(parallelPoints.pointPoints(), i)
226  {
227  if (master[i] != -1)
228  {
229  // I am master or slave
230  sharedPointLabels[nMaster] = i;
231  sharedPointAddr[nMaster] = master[i];
232  nMaster++;
233  }
234  }
235 
236  if (debug)
237  {
238  Pout<< "globalMeshData : nGlobalPoints_:" << nGlobalPoints_ << nl
239  << "globalMeshData : sharedPointLabels_:"
240  << sharedPointLabelsPtr_().size() << nl
241  << "globalMeshData : sharedPointAddr_:"
242  << sharedPointAddrPtr_().size() << endl;
243  }
244 }
245 
246 
247 void Foam::globalMeshData::countSharedEdges
248 (
249  const EdgeMap<labelList>& procSharedEdges,
250  EdgeMap<label>& globalShared,
251  label& sharedEdgeI
252 )
253 {
254  // Count occurrences of procSharedEdges in global shared edges table.
255  forAllConstIter(EdgeMap<labelList>, procSharedEdges, iter)
256  {
257  const edge& e = iter.key();
258 
259  EdgeMap<label>::iterator globalFnd = globalShared.find(e);
260 
261  if (globalFnd == globalShared.end())
262  {
263  // First time occurrence of this edge. Check how many we are adding.
264  if (iter().size() == 1)
265  {
266  // Only one edge. Mark with special value.
267  globalShared.insert(e, -1);
268  }
269  else
270  {
271  // Edge used more than once (even by local shared edges alone)
272  // so allocate proper shared edge label.
273  globalShared.insert(e, sharedEdgeI++);
274  }
275  }
276  else
277  {
278  if (globalFnd() == -1)
279  {
280  // Second time occurrence of this edge. Assign proper
281  // edge label.
282  globalFnd() = sharedEdgeI++;
283  }
284  }
285  }
286 }
287 
288 
289 void Foam::globalMeshData::calcSharedEdges() const
290 {
291  // Shared edges are shared between multiple processors. By their nature both
292  // of their endpoints are shared points. (but not all edges using two shared
293  // points are shared edges! There might e.g. be an edge between two
294  // unrelated clusters of shared points)
295 
296  if
297  (
298  nGlobalEdges_ != -1
299  || sharedEdgeLabelsPtr_.valid()
300  || sharedEdgeAddrPtr_.valid()
301  )
302  {
304  << "Shared edge addressing already done" << abort(FatalError);
305  }
306 
307 
308  const labelList& sharedPtAddr = sharedPointAddr();
309  const labelList& sharedPtLabels = sharedPointLabels();
310 
311  // Since don't want to construct pointEdges for whole mesh create
312  // Map for all shared points.
313  Map<label> meshToShared(2*sharedPtLabels.size());
314  forAll(sharedPtLabels, i)
315  {
316  meshToShared.insert(sharedPtLabels[i], i);
317  }
318 
319  // Find edges using shared points. Store correspondence to local edge
320  // numbering. Note that multiple local edges can have the same shared
321  // points! (for cyclics or separated processor patches)
322  EdgeMap<labelList> localShared(2*sharedPtAddr.size());
323 
324  const edgeList& edges = mesh_.edges();
325 
326  forAll(edges, edgeI)
327  {
328  const edge& e = edges[edgeI];
329 
330  Map<label>::const_iterator e0Fnd = meshToShared.find(e[0]);
331 
332  if (e0Fnd != meshToShared.end())
333  {
334  Map<label>::const_iterator e1Fnd = meshToShared.find(e[1]);
335 
336  if (e1Fnd != meshToShared.end())
337  {
338  // Found edge which uses shared points. Probably shared.
339 
340  // Construct the edge in shared points (or rather global indices
341  // of the shared points)
342  edge sharedEdge
343  (
344  sharedPtAddr[e0Fnd()],
345  sharedPtAddr[e1Fnd()]
346  );
347 
349  localShared.find(sharedEdge);
350 
351  if (iter == localShared.end())
352  {
353  // First occurrence of this point combination. Store.
354  localShared.insert(sharedEdge, labelList(1, edgeI));
355  }
356  else
357  {
358  // Add this edge to list of edge labels.
359  labelList& edgeLabels = iter();
360 
361  label sz = edgeLabels.size();
362  edgeLabels.setSize(sz+1);
363  edgeLabels[sz] = edgeI;
364  }
365  }
366  }
367  }
368 
369 
370  // Now we have a table on every processors which gives its edges which use
371  // shared points. Send this all to the master and have it allocate
372  // global edge numbers for it. But only allocate a global edge number for
373  // edge if it is used more than once!
374  // Note that we are now sending the whole localShared to the master whereas
375  // we only need the local count (i.e. the number of times a global edge is
376  // used). But then this only gets done once so not too bothered about the
377  // extra global communication.
378 
379  EdgeMap<label> globalShared(nGlobalPoints());
380 
381  if (Pstream::master())
382  {
383  label sharedEdgeI = 0;
384 
385  // Merge my shared edges into the global list
386  if (debug)
387  {
388  Pout<< "globalMeshData::calcSharedEdges : Merging in from proc0 : "
389  << localShared.size() << endl;
390  }
391  countSharedEdges(localShared, globalShared, sharedEdgeI);
392 
393  // Receive data from slaves and insert
394  if (Pstream::parRun())
395  {
396  for
397  (
398  int slave=Pstream::firstSlave();
399  slave<=Pstream::lastSlave();
400  slave++
401  )
402  {
403  // Receive the edges using shared points from the slave.
404  IPstream fromSlave(Pstream::commsTypes::blocking, slave);
405  EdgeMap<labelList> procSharedEdges(fromSlave);
406 
407  if (debug)
408  {
409  Pout<< "globalMeshData::calcSharedEdges : "
410  << "Merging in from proc"
411  << Foam::name(slave) << " : " << procSharedEdges.size()
412  << endl;
413  }
414  countSharedEdges(procSharedEdges, globalShared, sharedEdgeI);
415  }
416  }
417 
418  // Now our globalShared should have some edges with -1 as edge label
419  // These were only used once so are not proper shared edges.
420  // Remove them.
421  {
422  EdgeMap<label> oldSharedEdges(globalShared);
423 
424  globalShared.clear();
425 
426  forAllConstIter(EdgeMap<label>, oldSharedEdges, iter)
427  {
428  if (iter() != -1)
429  {
430  globalShared.insert(iter.key(), iter());
431  }
432  }
433  if (debug)
434  {
435  Pout<< "globalMeshData::calcSharedEdges : Filtered "
436  << oldSharedEdges.size()
437  << " down to " << globalShared.size() << endl;
438  }
439  }
440 
441 
442  // Send back to slaves.
443  if (Pstream::parRun())
444  {
445  for
446  (
447  int slave=Pstream::firstSlave();
448  slave<=Pstream::lastSlave();
449  slave++
450  )
451  {
452  // Receive the edges using shared points from the slave.
453  OPstream toSlave(Pstream::commsTypes::blocking, slave);
454  toSlave << globalShared;
455  }
456  }
457  }
458  else
459  {
460  // Send local edges to master
461  {
462  OPstream toMaster
463  (
466  );
467  toMaster << localShared;
468  }
469  // Receive merged edges from master.
470  {
471  IPstream fromMaster
472  (
475  );
476  fromMaster >> globalShared;
477  }
478  }
479 
480  // Now use the global shared edges list (globalShared) to classify my local
481  // ones (localShared)
482 
483  nGlobalEdges_ = globalShared.size();
484 
485  DynamicList<label> dynSharedEdgeLabels(globalShared.size());
486  DynamicList<label> dynSharedEdgeAddr(globalShared.size());
487 
488  forAllConstIter(EdgeMap<labelList>, localShared, iter)
489  {
490  const edge& e = iter.key();
491 
492  EdgeMap<label>::const_iterator edgeFnd = globalShared.find(e);
493 
494  if (edgeFnd != globalShared.end())
495  {
496  // My local edge is indeed a shared one. Go through all local edge
497  // labels with this point combination.
498  const labelList& edgeLabels = iter();
499 
500  forAll(edgeLabels, i)
501  {
502  // Store label of local mesh edge
503  dynSharedEdgeLabels.append(edgeLabels[i]);
504 
505  // Store label of shared edge
506  dynSharedEdgeAddr.append(edgeFnd());
507  }
508  }
509  }
510 
511  sharedEdgeLabelsPtr_.reset(new labelList());
512  labelList& sharedEdgeLabels = sharedEdgeLabelsPtr_();
513  sharedEdgeLabels.transfer(dynSharedEdgeLabels);
514 
515  sharedEdgeAddrPtr_.reset(new labelList());
516  labelList& sharedEdgeAddr = sharedEdgeAddrPtr_();
517  sharedEdgeAddr.transfer(dynSharedEdgeAddr);
518 
519  if (debug)
520  {
521  Pout<< "globalMeshData : nGlobalEdges_:" << nGlobalEdges_ << nl
522  << "globalMeshData : sharedEdgeLabels:" << sharedEdgeLabels.size()
523  << nl
524  << "globalMeshData : sharedEdgeAddr:" << sharedEdgeAddr.size()
525  << endl;
526  }
527 }
528 
529 
530 void Foam::globalMeshData::calcGlobalPointSlaves() const
531 {
532  if (debug)
533  {
534  Pout<< "globalMeshData::calcGlobalPointSlaves() :"
535  << " calculating coupled master to slave point addressing."
536  << endl;
537  }
538 
539  // Calculate connected points for master points.
540  globalPoints globalData(mesh_, coupledPatch(), true, true);
541 
542  globalPointSlavesPtr_.reset
543  (
544  new labelListList
545  (
546  move(globalData.pointPoints())
547  )
548  );
549  globalPointTransformedSlavesPtr_.reset
550  (
551  new labelListList
552  (
553  move(globalData.transformedPointPoints())
554  )
555  );
556 
557  globalPointSlavesMapPtr_.reset
558  (
559  new mapDistribute
560  (
561  move(globalData.map())
562  )
563  );
564 }
565 
566 
567 void Foam::globalMeshData::calcPointConnectivity
568 (
569  List<labelPairList>& allPointConnectivity
570 ) const
571 {
572  const globalIndexAndTransform& transforms = globalTransforms();
573  const labelListList& slaves = globalPointSlaves();
574  const labelListList& transformedSlaves = globalPointTransformedSlaves();
575 
576 
577  // Create field with my local data
578  labelPairList myData(globalPointSlavesMap().constructSize());
579  forAll(slaves, pointi)
580  {
581  myData[pointi] = transforms.encode
582  (
584  pointi,
585  transforms.nullTransformIndex()
586  );
587  }
588  // Send to master
589  globalPointSlavesMap().distribute(myData);
590 
591 
592  // String of connected points with their transform
593  allPointConnectivity.setSize(globalPointSlavesMap().constructSize());
594  allPointConnectivity = labelPairList(0);
595 
596  // Pass1: do the master points since these also update local slaves
597  // (e.g. from local cyclics)
598  forAll(slaves, pointi)
599  {
600  // Reconstruct string of connected points
601  const labelList& pSlaves = slaves[pointi];
602  const labelList& pTransformSlaves = transformedSlaves[pointi];
603 
604  if (pSlaves.size()+pTransformSlaves.size())
605  {
606  labelPairList& pConnectivity = allPointConnectivity[pointi];
607 
608  pConnectivity.setSize(1+pSlaves.size()+pTransformSlaves.size());
609  label connI = 0;
610 
611  // Add myself
612  pConnectivity[connI++] = myData[pointi];
613  // Add untransformed points
614  forAll(pSlaves, i)
615  {
616  pConnectivity[connI++] = myData[pSlaves[i]];
617  }
618  // Add transformed points.
619  forAll(pTransformSlaves, i)
620  {
621  // Get transform from index
622  label transformI = globalPointSlavesMap().whichTransform
623  (
624  pTransformSlaves[i]
625  );
626  // Add transform to connectivity
627  const labelPair& n = myData[pTransformSlaves[i]];
628  label proci = transforms.processor(n);
629  label index = transforms.index(n);
630  pConnectivity[connI++] = transforms.encode
631  (
632  proci,
633  index,
634  transformI
635  );
636  }
637 
638  // Put back in slots
639  forAll(pSlaves, i)
640  {
641  allPointConnectivity[pSlaves[i]] = pConnectivity;
642  }
643  forAll(pTransformSlaves, i)
644  {
645  allPointConnectivity[pTransformSlaves[i]] = pConnectivity;
646  }
647  }
648  }
649 
650 
651  // Pass2: see if anything is still unset (should not be the case)
652  forAll(slaves, pointi)
653  {
654  labelPairList& pConnectivity = allPointConnectivity[pointi];
655 
656  if (pConnectivity.size() == 0)
657  {
658  pConnectivity.setSize(1, myData[pointi]);
659  }
660  }
661 
662 
663  globalPointSlavesMap().reverseDistribute
664  (
665  slaves.size(),
666  allPointConnectivity
667  );
668 }
669 
670 
671 void Foam::globalMeshData::calcGlobalPointEdges
672 (
673  labelListList& globalPointEdges,
674  List<labelPairList>& globalPointPoints
675 ) const
676 {
677  const edgeList& edges = coupledPatch().edges();
678  const labelListList& pointEdges = coupledPatch().pointEdges();
679  const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
680  const labelListList& slaves = globalPointSlaves();
681  const labelListList& transformedSlaves = globalPointTransformedSlaves();
682  const globalIndexAndTransform& transforms = globalTransforms();
683 
684 
685  // Create local version
686  globalPointEdges.setSize(globalPointSlavesMap().constructSize());
687  globalPointPoints.setSize(globalPointSlavesMap().constructSize());
688  forAll(pointEdges, pointi)
689  {
690  const labelList& pEdges = pointEdges[pointi];
691  labelList& globalPEdges = globalPointEdges[pointi];
692  globalPEdges.setSize(pEdges.size());
693  forAll(pEdges, i)
694  {
695  globalPEdges[i] = globalEdgeNumbers.toGlobal(pEdges[i]);
696  }
697 
698  labelPairList& globalPPoints = globalPointPoints[pointi];
699  globalPPoints.setSize(pEdges.size());
700  forAll(pEdges, i)
701  {
702  label otherPointi = edges[pEdges[i]].otherVertex(pointi);
703  globalPPoints[i] = transforms.encode
704  (
706  otherPointi,
707  transforms.nullTransformIndex()
708  );
709  }
710  }
711 
712  // Pull slave data to master. Dummy transform.
713  globalPointSlavesMap().distribute(globalPointEdges);
714  globalPointSlavesMap().distribute(globalPointPoints);
715  // Add all pointEdges
716  forAll(slaves, pointi)
717  {
718  const labelList& pSlaves = slaves[pointi];
719  const labelList& pTransformSlaves = transformedSlaves[pointi];
720 
721  label n = 0;
722  forAll(pSlaves, i)
723  {
724  n += globalPointEdges[pSlaves[i]].size();
725  }
726  forAll(pTransformSlaves, i)
727  {
728  n += globalPointEdges[pTransformSlaves[i]].size();
729  }
730 
731  // Add all the point edges of the slaves to those of the (master) point
732  {
733  labelList& globalPEdges = globalPointEdges[pointi];
734  label sz = globalPEdges.size();
735  globalPEdges.setSize(sz+n);
736  forAll(pSlaves, i)
737  {
738  const labelList& otherData = globalPointEdges[pSlaves[i]];
739  forAll(otherData, j)
740  {
741  globalPEdges[sz++] = otherData[j];
742  }
743  }
744  forAll(pTransformSlaves, i)
745  {
746  const labelList& otherData =
747  globalPointEdges[pTransformSlaves[i]];
748  forAll(otherData, j)
749  {
750  globalPEdges[sz++] = otherData[j];
751  }
752  }
753 
754  // Put back in slots
755  forAll(pSlaves, i)
756  {
757  globalPointEdges[pSlaves[i]] = globalPEdges;
758  }
759  forAll(pTransformSlaves, i)
760  {
761  globalPointEdges[pTransformSlaves[i]] = globalPEdges;
762  }
763  }
764 
765 
766  // Same for corresponding pointPoints
767  {
768  labelPairList& globalPPoints = globalPointPoints[pointi];
769  label sz = globalPPoints.size();
770  globalPPoints.setSize(sz + n);
771 
772  // Add untransformed points
773  forAll(pSlaves, i)
774  {
775  const labelPairList& otherData = globalPointPoints[pSlaves[i]];
776  forAll(otherData, j)
777  {
778  globalPPoints[sz++] = otherData[j];
779  }
780  }
781  // Add transformed points.
782  forAll(pTransformSlaves, i)
783  {
784  // Get transform from index
785  label transformI = globalPointSlavesMap().whichTransform
786  (
787  pTransformSlaves[i]
788  );
789 
790  const labelPairList& otherData =
791  globalPointPoints[pTransformSlaves[i]];
792  forAll(otherData, j)
793  {
794  // Add transform to connectivity
795  const labelPair& n = otherData[j];
796  label proci = transforms.processor(n);
797  label index = transforms.index(n);
798  globalPPoints[sz++] = transforms.encode
799  (
800  proci,
801  index,
802  transformI
803  );
804  }
805  }
806 
807  // Put back in slots
808  forAll(pSlaves, i)
809  {
810  globalPointPoints[pSlaves[i]] = globalPPoints;
811  }
812  forAll(pTransformSlaves, i)
813  {
814  globalPointPoints[pTransformSlaves[i]] = globalPPoints;
815  }
816  }
817  }
818  // Push back
819  globalPointSlavesMap().reverseDistribute
820  (
821  slaves.size(),
822  globalPointEdges
823  );
824  // Push back
825  globalPointSlavesMap().reverseDistribute
826  (
827  slaves.size(),
828  globalPointPoints
829  );
830 }
831 
832 
833 Foam::label Foam::globalMeshData::findTransform
834 (
835  const labelPairList& info,
836  const labelPair& remotePoint,
837  const label localPoint
838 ) const
839 {
840  const globalIndexAndTransform& transforms = globalTransforms();
841 
842  const label remoteProci = transforms.processor(remotePoint);
843  const label remoteIndex = transforms.index(remotePoint);
844 
845  label remoteTransformI = -1;
846  label localTransformI = -1;
847  forAll(info, i)
848  {
849  label proci = transforms.processor(info[i]);
850  label pointi = transforms.index(info[i]);
851  label transformI = transforms.transformIndex(info[i]);
852 
853  if (proci == Pstream::myProcNo() && pointi == localPoint)
854  {
855  localTransformI = transformI;
856  // Pout<< "For local :" << localPoint
857  // << " found transform:" << localTransformI
858  // << endl;
859  }
860  if (proci == remoteProci && pointi == remoteIndex)
861  {
862  remoteTransformI = transformI;
863  // Pout<< "For remote:" << remotePoint
864  // << " found transform:" << remoteTransformI
865  // << " at index:" << i
866  // << endl;
867  }
868  }
869 
870  if (remoteTransformI == -1 || localTransformI == -1)
871  {
873  << "Problem. Cannot find " << remotePoint
874  << " or " << localPoint << " "
875  << coupledPatch().localPoints()[localPoint]
876  << " in " << info
877  << endl
878  << "remoteTransformI:" << remoteTransformI << endl
879  << "localTransformI:" << localTransformI
880  << abort(FatalError);
881  }
882 
883  return transforms.subtractTransformIndex
884  (
885  remoteTransformI,
886  localTransformI
887  );
888 }
889 
890 
891 void Foam::globalMeshData::calcGlobalEdgeSlaves() const
892 {
893  if (debug)
894  {
895  Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
896  << " calculating coupled master to slave edge addressing." << endl;
897  }
898 
899  const edgeList& edges = coupledPatch().edges();
900  const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
901  const globalIndexAndTransform& transforms = globalTransforms();
902 
903 
904  // The whole problem with deducting edge-connectivity from
905  // point-connectivity is that one of the endpoints might be
906  // a local master but the other endpoint might not. So we first
907  // need to make sure that all points know about connectivity and
908  // the transformations.
909 
910 
911  // 1. collect point connectivity - basically recreating globalPoints output.
912  // All points will now have a string of coupled points. The transforms are
913  // in respect to the master.
914  List<labelPairList> allPointConnectivity;
915  calcPointConnectivity(allPointConnectivity);
916 
917 
918  // 2. Get all pointEdges and pointPoints
919  // Coupled point to global coupled edges and corresponding endpoint.
920  labelListList globalPointEdges;
921  List<labelPairList> globalPointPoints;
922  calcGlobalPointEdges(globalPointEdges, globalPointPoints);
923 
924 
925  // 3. Now all points have
926  // - all the connected points with original transform
927  // - all the connected global edges
928 
929  // Now all we need to do is go through all the edges and check
930  // both endpoints. If there is a edge between the two which is
931  // produced by transforming both points in the same way it is a shared
932  // edge.
933 
934  // Collect strings of connected edges.
935  List<labelPairList> allEdgeConnectivity(edges.size());
936 
937  forAll(edges, edgeI)
938  {
939  const edge& e = edges[edgeI];
940  const labelList& pEdges0 = globalPointEdges[e[0]];
941  const labelPairList& pPoints0 = globalPointPoints[e[0]];
942  const labelList& pEdges1 = globalPointEdges[e[1]];
943  const labelPairList& pPoints1 = globalPointPoints[e[1]];
944 
945  // Most edges will be size 2
946  DynamicList<labelPair> eEdges(2);
947  // Append myself.
948  eEdges.append
949  (
950  transforms.encode
951  (
953  edgeI,
954  transforms.nullTransformIndex()
955  )
956  );
957 
958  forAll(pEdges0, i)
959  {
960  forAll(pEdges1, j)
961  {
962  if
963  (
964  pEdges0[i] == pEdges1[j]
965  && pEdges0[i] != globalEdgeNumbers.toGlobal(edgeI)
966  )
967  {
968  // Found a shared edge. Now check if the endpoints
969  // go through the same transformation.
970  // Local: e[0] remote:pPoints1[j]
971  // Local: e[1] remote:pPoints0[i]
972 
973 
974  // Find difference in transforms to go from point on remote
975  // edge (pPoints1[j]) to this point.
976 
977  label transform0 = findTransform
978  (
979  allPointConnectivity[e[0]],
980  pPoints1[j],
981  e[0]
982  );
983  label transform1 = findTransform
984  (
985  allPointConnectivity[e[1]],
986  pPoints0[i],
987  e[1]
988  );
989 
990  if (transform0 == transform1)
991  {
992  label proci = globalEdgeNumbers.whichProcID(pEdges0[i]);
993  eEdges.append
994  (
995  transforms.encode
996  (
997  proci,
998  globalEdgeNumbers.toLocal(proci, pEdges0[i]),
999  transform0
1000  )
1001  );
1002  }
1003  }
1004  }
1005  }
1006 
1007  allEdgeConnectivity[edgeI].transfer(eEdges);
1008  sort
1009  (
1010  allEdgeConnectivity[edgeI],
1011  globalIndexAndTransform::less(transforms)
1012  );
1013  }
1014 
1015  // We now have - in allEdgeConnectivity - a list of edges which are shared
1016  // between multiple processors. Filter into non-transformed and transformed
1017  // connections.
1018 
1019  globalEdgeSlavesPtr_.reset(new labelListList(edges.size()));
1020  labelListList& globalEdgeSlaves = globalEdgeSlavesPtr_();
1021  List<labelPairList> transformedEdges(edges.size());
1022  forAll(allEdgeConnectivity, edgeI)
1023  {
1024  const labelPairList& edgeInfo = allEdgeConnectivity[edgeI];
1025  if (edgeInfo.size() >= 2)
1026  {
1027  const labelPair& masterInfo = edgeInfo[0];
1028 
1029  // Check if master edge (= first element (since sorted)) is me.
1030  if
1031  (
1032  (
1033  transforms.processor(masterInfo)
1034  == Pstream::myProcNo()
1035  )
1036  && (transforms.index(masterInfo) == edgeI)
1037  )
1038  {
1039  // Sort into transformed and untransformed
1040  labelList& eEdges = globalEdgeSlaves[edgeI];
1041  eEdges.setSize(edgeInfo.size()-1);
1042 
1043  labelPairList& trafoEEdges = transformedEdges[edgeI];
1044  trafoEEdges.setSize(edgeInfo.size()-1);
1045 
1046  label nonTransformI = 0;
1047  label transformI = 0;
1048 
1049  for (label i = 1; i < edgeInfo.size(); i++)
1050  {
1051  const labelPair& info = edgeInfo[i];
1052  label proci = transforms.processor(info);
1053  label index = transforms.index(info);
1054  label transform = transforms.transformIndex
1055  (
1056  info
1057  );
1058 
1059  if (transform == transforms.nullTransformIndex())
1060  {
1061  eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
1062  (
1063  proci,
1064  index
1065  );
1066  }
1067  else
1068  {
1069  trafoEEdges[transformI++] = info;
1070  }
1071  }
1072 
1073  eEdges.setSize(nonTransformI);
1074  trafoEEdges.setSize(transformI);
1075  }
1076  }
1077  }
1078 
1079 
1080  // Construct map
1081  globalEdgeTransformedSlavesPtr_.reset(new labelListList());
1082 
1083  List<Map<label>> compactMap(Pstream::nProcs());
1084  globalEdgeSlavesMapPtr_.reset
1085  (
1086  new mapDistribute
1087  (
1088  globalEdgeNumbers,
1089  globalEdgeSlaves,
1090 
1091  transforms,
1092  transformedEdges,
1093  globalEdgeTransformedSlavesPtr_(),
1094 
1095  compactMap
1096  )
1097  );
1098 
1099 
1100  if (debug)
1101  {
1102  Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
1103  << " coupled edges:" << edges.size()
1104  << " additional coupled edges:"
1105  << globalEdgeSlavesMapPtr_().constructSize() - edges.size()
1106  << endl;
1107  }
1108 }
1109 
1110 
1111 void Foam::globalMeshData::calcGlobalEdgeOrientation() const
1112 {
1113  if (debug)
1114  {
1115  Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
1116  << " calculating edge orientation w.r.t. master edge." << endl;
1117  }
1118 
1119  const globalIndex& globalPoints = globalPointNumbering();
1120 
1121  // 1. Determine master point
1122  labelList masterPoint;
1123  {
1124  const mapDistribute& map = globalPointSlavesMap();
1125 
1126  masterPoint.setSize(map.constructSize());
1127  masterPoint = labelMax;
1128 
1129  for (label pointi = 0; pointi < coupledPatch().nPoints(); pointi++)
1130  {
1131  masterPoint[pointi] = globalPoints.toGlobal(pointi);
1132  }
1133  syncData
1134  (
1135  masterPoint,
1136  globalPointSlaves(),
1137  globalPointTransformedSlaves(),
1138  map,
1139  minEqOp<label>()
1140  );
1141  }
1142 
1143  // Now all points should know who is master by comparing their global
1144  // pointID with the masterPointID. We now can use this information
1145  // to find the orientation of the master edge.
1146 
1147  {
1148  const mapDistribute& map = globalEdgeSlavesMap();
1149  const labelListList& slaves = globalEdgeSlaves();
1150  const labelListList& transformedSlaves = globalEdgeTransformedSlaves();
1151 
1152  // Distribute orientation of master edge (in masterPoint numbering)
1153  labelPairList masterEdgeVerts(map.constructSize());
1154  masterEdgeVerts = labelPair(labelMax, labelMax);
1155 
1156  for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++)
1157  {
1158  if
1159  (
1160  (
1161  slaves[edgeI].size()
1162  + transformedSlaves[edgeI].size()
1163  )
1164  > 0
1165  )
1166  {
1167  // I am master. Fill in my masterPoint equivalent.
1168 
1169  const edge& e = coupledPatch().edges()[edgeI];
1170  masterEdgeVerts[edgeI] = labelPair
1171  (
1172  masterPoint[e[0]],
1173  masterPoint[e[1]]
1174  );
1175  }
1176  }
1177  syncData
1178  (
1179  masterEdgeVerts,
1180  slaves,
1181  transformedSlaves,
1182  map,
1184  );
1185 
1186  // Now check my edges on how they relate to the master's edgeVerts
1187  globalEdgeOrientationPtr_.reset
1188  (
1189  new PackedBoolList(coupledPatch().nEdges())
1190  );
1191  PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_();
1192 
1193  forAll(coupledPatch().edges(), edgeI)
1194  {
1195  const edge& e = coupledPatch().edges()[edgeI];
1196  const labelPair masterE
1197  (
1198  masterPoint[e[0]],
1199  masterPoint[e[1]]
1200  );
1201 
1202  label stat = labelPair::compare
1203  (
1204  masterE,
1205  masterEdgeVerts[edgeI]
1206  );
1207  if (stat == 0)
1208  {
1210  << "problem : my edge:" << e
1211  << " in master points:" << masterE
1212  << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI]
1213  << exit(FatalError);
1214  }
1215  else
1216  {
1217  globalEdgeOrientation[edgeI] = (stat == 1);
1218  }
1219  }
1220  }
1221 
1222  if (debug)
1223  {
1224  Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
1225  << " finished calculating edge orientation."
1226  << endl;
1227  }
1228 }
1229 
1230 
1231 void Foam::globalMeshData::calcPointBoundaryFaces
1232 (
1233  labelListList& pointBoundaryFaces
1234 ) const
1235 {
1236  const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
1237  const Map<label>& meshPointMap = coupledPatch().meshPointMap();
1238 
1239  // 1. Count
1240 
1241  labelList nPointFaces(coupledPatch().nPoints(), 0);
1242 
1243  forAll(bMesh, patchi)
1244  {
1245  const polyPatch& pp = bMesh[patchi];
1246 
1247  if (!pp.coupled())
1248  {
1249  forAll(pp, i)
1250  {
1251  const face& f = pp[i];
1252 
1253  forAll(f, fp)
1254  {
1255  Map<label>::const_iterator iter = meshPointMap.find
1256  (
1257  f[fp]
1258  );
1259  if (iter != meshPointMap.end())
1260  {
1261  nPointFaces[iter()]++;
1262  }
1263  }
1264  }
1265  }
1266  }
1267 
1268 
1269  // 2. Size
1270 
1271  pointBoundaryFaces.setSize(coupledPatch().nPoints());
1272  forAll(nPointFaces, pointi)
1273  {
1274  pointBoundaryFaces[pointi].setSize(nPointFaces[pointi]);
1275  }
1276  nPointFaces = 0;
1277 
1278 
1279  // 3. Fill
1280 
1281  forAll(bMesh, patchi)
1282  {
1283  const polyPatch& pp = bMesh[patchi];
1284 
1285  if (!pp.coupled())
1286  {
1287  forAll(pp, i)
1288  {
1289  const face& f = pp[i];
1290  forAll(f, fp)
1291  {
1292  Map<label>::const_iterator iter = meshPointMap.find
1293  (
1294  f[fp]
1295  );
1296  if (iter != meshPointMap.end())
1297  {
1298  label bFacei =
1299  pp.start() + i - mesh_.nInternalFaces();
1300  pointBoundaryFaces[iter()][nPointFaces[iter()]++] =
1301  bFacei;
1302  }
1303  }
1304  }
1305  }
1306  }
1307 }
1308 
1309 
1310 void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
1311 {
1312  if (debug)
1313  {
1314  Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
1315  << " calculating coupled point to boundary face addressing."
1316  << endl;
1317  }
1318 
1319  // Construct local point to (uncoupled)boundaryfaces.
1320  labelListList pointBoundaryFaces;
1321  calcPointBoundaryFaces(pointBoundaryFaces);
1322 
1323 
1324  // Global indices for boundary faces
1325  globalBoundaryFaceNumberingPtr_.reset
1326  (
1327  new globalIndex(mesh_.nFaces()-mesh_.nInternalFaces())
1328  );
1329  globalIndex& globalIndices = globalBoundaryFaceNumberingPtr_();
1330 
1331 
1332  // Convert local boundary faces to global numbering
1333  globalPointBoundaryFacesPtr_.reset
1334  (
1335  new labelListList(globalPointSlavesMap().constructSize())
1336  );
1337  labelListList& globalPointBoundaryFaces = globalPointBoundaryFacesPtr_();
1338 
1339  forAll(pointBoundaryFaces, pointi)
1340  {
1341  const labelList& bFaces = pointBoundaryFaces[pointi];
1342  labelList& globalFaces = globalPointBoundaryFaces[pointi];
1343  globalFaces.setSize(bFaces.size());
1344  forAll(bFaces, i)
1345  {
1346  globalFaces[i] = globalIndices.toGlobal(bFaces[i]);
1347  }
1348  }
1349 
1350 
1351  // Pull slave pointBoundaryFaces to master
1352  globalPointSlavesMap().distribute
1353  (
1354  globalPointBoundaryFaces,
1355  true // put data on transformed points into correct slots
1356  );
1357 
1358 
1359  // Merge slave labels into master globalPointBoundaryFaces.
1360  // Split into untransformed and transformed values.
1361  const labelListList& pointSlaves = globalPointSlaves();
1362  const labelListList& pointTransformSlaves =
1363  globalPointTransformedSlaves();
1364  const globalIndexAndTransform& transforms = globalTransforms();
1365 
1366 
1367  // Any faces coming in through transformation
1368  List<labelPairList> transformedFaces(pointSlaves.size());
1369 
1370 
1371  forAll(pointSlaves, pointi)
1372  {
1373  const labelList& slaves = pointSlaves[pointi];
1374  const labelList& transformedSlaves = pointTransformSlaves[pointi];
1375 
1376  if (slaves.size() > 0)
1377  {
1378  labelList& myBFaces = globalPointBoundaryFaces[pointi];
1379  label sz = myBFaces.size();
1380 
1381  // Count
1382  label n = 0;
1383  forAll(slaves, i)
1384  {
1385  n += globalPointBoundaryFaces[slaves[i]].size();
1386  }
1387  // Fill
1388  myBFaces.setSize(sz+n);
1389  n = sz;
1390  forAll(slaves, i)
1391  {
1392  const labelList& slaveBFaces =
1393  globalPointBoundaryFaces[slaves[i]];
1394 
1395  // Add all slaveBFaces. Note that need to check for
1396  // uniqueness only in case of cyclics.
1397 
1398  forAll(slaveBFaces, j)
1399  {
1400  label slave = slaveBFaces[j];
1401  if (findIndex(SubList<label>(myBFaces, sz), slave) == -1)
1402  {
1403  myBFaces[n++] = slave;
1404  }
1405  }
1406  }
1407  myBFaces.setSize(n);
1408  }
1409 
1410 
1411  if (transformedSlaves.size() > 0)
1412  {
1413  const labelList& untrafoFaces = globalPointBoundaryFaces[pointi];
1414 
1415  labelPairList& myBFaces = transformedFaces[pointi];
1416  label sz = myBFaces.size();
1417 
1418  // Count
1419  label n = 0;
1420  forAll(transformedSlaves, i)
1421  {
1422  n += globalPointBoundaryFaces[transformedSlaves[i]].size();
1423  }
1424  // Fill
1425  myBFaces.setSize(sz+n);
1426  n = sz;
1427  forAll(transformedSlaves, i)
1428  {
1429  label transformI = globalPointSlavesMap().whichTransform
1430  (
1431  transformedSlaves[i]
1432  );
1433 
1434  const labelList& slaveBFaces =
1435  globalPointBoundaryFaces[transformedSlaves[i]];
1436 
1437  forAll(slaveBFaces, j)
1438  {
1439  label slave = slaveBFaces[j];
1440  // Check that same face not already present untransformed
1441  if (findIndex(untrafoFaces, slave)== -1)
1442  {
1443  label proci = globalIndices.whichProcID(slave);
1444  label facei = globalIndices.toLocal(proci, slave);
1445 
1446  myBFaces[n++] = transforms.encode
1447  (
1448  proci,
1449  facei,
1450  transformI
1451  );
1452  }
1453  }
1454  }
1455  myBFaces.setSize(n);
1456  }
1457 
1458 
1459  if (slaves.size() + transformedSlaves.size() == 0)
1460  {
1461  globalPointBoundaryFaces[pointi].clear();
1462  }
1463  }
1464 
1465  // Construct a map to get the face data directly
1466  List<Map<label>> compactMap(Pstream::nProcs());
1467 
1468  globalPointTransformedBoundaryFacesPtr_.reset
1469  (
1470  new labelListList(transformedFaces.size())
1471  );
1472 
1473  globalPointBoundaryFacesMapPtr_.reset
1474  (
1475  new mapDistribute
1476  (
1477  globalIndices,
1478  globalPointBoundaryFaces,
1479 
1480  transforms,
1481  transformedFaces,
1482  globalPointTransformedBoundaryFacesPtr_(),
1483 
1484  compactMap
1485  )
1486  );
1487  globalPointBoundaryFaces.setSize(coupledPatch().nPoints());
1488  globalPointTransformedBoundaryFacesPtr_().setSize(coupledPatch().nPoints());
1489 
1490  if (debug)
1491  {
1492  Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
1493  << " coupled points:" << coupledPatch().nPoints()
1494  << " local boundary faces:" << globalIndices.localSize()
1495  << " additional coupled faces:"
1496  << globalPointBoundaryFacesMapPtr_().constructSize()
1497  - globalIndices.localSize()
1498  << endl;
1499  }
1500 }
1501 
1502 
1503 void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
1504 {
1505  if (debug)
1506  {
1507  Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
1508  << " calculating coupled point to boundary cell addressing."
1509  << endl;
1510  }
1511 
1512  // Create map of boundary cells and point-cell addressing
1513  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1514 
1515  label bCelli = 0;
1516  Map<label> meshCellMap(4*coupledPatch().nPoints());
1517  DynamicList<label> cellMap(meshCellMap.size());
1518 
1519  // Create addressing for point to boundary cells (local)
1520  labelListList pointBoundaryCells(coupledPatch().nPoints());
1521 
1522  forAll(coupledPatch().meshPoints(), pointi)
1523  {
1524  label meshPointi = coupledPatch().meshPoints()[pointi];
1525  const labelList& pCells = mesh_.pointCells(meshPointi);
1526 
1527  labelList& bCells = pointBoundaryCells[pointi];
1528  bCells.setSize(pCells.size());
1529 
1530  forAll(pCells, i)
1531  {
1532  label celli = pCells[i];
1533  Map<label>::iterator fnd = meshCellMap.find(celli);
1534 
1535  if (fnd != meshCellMap.end())
1536  {
1537  bCells[i] = fnd();
1538  }
1539  else
1540  {
1541  meshCellMap.insert(celli, bCelli);
1542  cellMap.append(celli);
1543  bCells[i] = bCelli;
1544  bCelli++;
1545  }
1546  }
1547  }
1548 
1549 
1550  boundaryCellsPtr_.reset(new labelList());
1551  labelList& boundaryCells = boundaryCellsPtr_();
1552  boundaryCells.transfer(cellMap.shrink());
1553 
1554 
1555  // Convert point-cells to global (boundary)cell numbers
1556  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1557 
1558  globalBoundaryCellNumberingPtr_.reset
1559  (
1560  new globalIndex(boundaryCells.size())
1561  );
1562  globalIndex& globalIndices = globalBoundaryCellNumberingPtr_();
1563 
1564 
1565  globalPointBoundaryCellsPtr_.reset
1566  (
1567  new labelListList(globalPointSlavesMap().constructSize())
1568  );
1569  labelListList& globalPointBoundaryCells = globalPointBoundaryCellsPtr_();
1570 
1571  forAll(pointBoundaryCells, pointi)
1572  {
1573  const labelList& pCells = pointBoundaryCells[pointi];
1574  labelList& globalCells = globalPointBoundaryCells[pointi];
1575  globalCells.setSize(pCells.size());
1576  forAll(pCells, i)
1577  {
1578  globalCells[i] = globalIndices.toGlobal(pCells[i]);
1579  }
1580  }
1581 
1582 
1583  // Pull slave pointBoundaryCells to master
1584  globalPointSlavesMap().distribute
1585  (
1586  globalPointBoundaryCells,
1587  true // put data on transformed points into correct slots
1588  );
1589 
1590 
1591  // Merge slave labels into master globalPointBoundaryCells
1592  const labelListList& pointSlaves = globalPointSlaves();
1593  const labelListList& pointTransformSlaves =
1594  globalPointTransformedSlaves();
1595  const globalIndexAndTransform& transforms = globalTransforms();
1596 
1597  List<labelPairList> transformedCells(pointSlaves.size());
1598 
1599 
1600  forAll(pointSlaves, pointi)
1601  {
1602  const labelList& slaves = pointSlaves[pointi];
1603  const labelList& transformedSlaves = pointTransformSlaves[pointi];
1604 
1605  if (slaves.size() > 0)
1606  {
1607  labelList& myBCells = globalPointBoundaryCells[pointi];
1608  label sz = myBCells.size();
1609 
1610  // Count
1611  label n = 0;
1612  forAll(slaves, i)
1613  {
1614  n += globalPointBoundaryCells[slaves[i]].size();
1615  }
1616  // Fill
1617  myBCells.setSize(sz+n);
1618  n = sz;
1619  forAll(slaves, i)
1620  {
1621  const labelList& slaveBCells =
1622  globalPointBoundaryCells[slaves[i]];
1623 
1624  // Add all slaveBCells. Note that need to check for
1625  // uniqueness only in case of cyclics.
1626 
1627  forAll(slaveBCells, j)
1628  {
1629  label slave = slaveBCells[j];
1630  if (findIndex(SubList<label>(myBCells, sz), slave) == -1)
1631  {
1632  myBCells[n++] = slave;
1633  }
1634  }
1635  }
1636  myBCells.setSize(n);
1637  }
1638 
1639 
1640  if (transformedSlaves.size() > 0)
1641  {
1642  const labelList& untrafoCells = globalPointBoundaryCells[pointi];
1643 
1644  labelPairList& myBCells = transformedCells[pointi];
1645  label sz = myBCells.size();
1646 
1647  // Count
1648  label n = 0;
1649  forAll(transformedSlaves, i)
1650  {
1651  n += globalPointBoundaryCells[transformedSlaves[i]].size();
1652  }
1653  // Fill
1654  myBCells.setSize(sz+n);
1655  n = sz;
1656  forAll(transformedSlaves, i)
1657  {
1658  label transformI = globalPointSlavesMap().whichTransform
1659  (
1660  transformedSlaves[i]
1661  );
1662 
1663  const labelList& slaveBCells =
1664  globalPointBoundaryCells[transformedSlaves[i]];
1665 
1666  forAll(slaveBCells, j)
1667  {
1668  label slave = slaveBCells[j];
1669 
1670  // Check that same cell not already present untransformed
1671  if (findIndex(untrafoCells, slave)== -1)
1672  {
1673  label proci = globalIndices.whichProcID(slave);
1674  label celli = globalIndices.toLocal(proci, slave);
1675  myBCells[n++] = transforms.encode
1676  (
1677  proci,
1678  celli,
1679  transformI
1680  );
1681  }
1682  }
1683  }
1684  myBCells.setSize(n);
1685  }
1686 
1687  if (slaves.size() + transformedSlaves.size() == 0)
1688  {
1689  globalPointBoundaryCells[pointi].clear();
1690  }
1691  }
1692 
1693  // Construct a map to get the cell data directly
1694  List<Map<label>> compactMap(Pstream::nProcs());
1695 
1696  globalPointTransformedBoundaryCellsPtr_.reset
1697  (
1698  new labelListList(transformedCells.size())
1699  );
1700 
1701  globalPointBoundaryCellsMapPtr_.reset
1702  (
1703  new mapDistribute
1704  (
1705  globalIndices,
1706  globalPointBoundaryCells,
1707 
1708  transforms,
1709  transformedCells,
1710  globalPointTransformedBoundaryCellsPtr_(),
1711 
1712  compactMap
1713  )
1714  );
1715  globalPointBoundaryCells.setSize(coupledPatch().nPoints());
1716  globalPointTransformedBoundaryCellsPtr_().setSize(coupledPatch().nPoints());
1717 
1718  if (debug)
1719  {
1720  Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
1721  << " coupled points:" << coupledPatch().nPoints()
1722  << " local boundary cells:" << globalIndices.localSize()
1723  << " additional coupled cells:"
1724  << globalPointBoundaryCellsMapPtr_().constructSize()
1725  - globalIndices.localSize()
1726  << endl;
1727  }
1728 }
1729 
1730 
1731 void Foam::globalMeshData::calcGlobalCoPointSlaves() const
1732 {
1733  if (debug)
1734  {
1735  Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
1736  << " calculating coupled master to collocated"
1737  << " slave point addressing." << endl;
1738  }
1739 
1740  // Calculate connected points for master points.
1741  globalPoints globalData(mesh_, coupledPatch(), true, false);
1742 
1743  globalCoPointSlavesPtr_.reset
1744  (
1745  new labelListList
1746  (
1747  move(globalData.pointPoints())
1748  )
1749  );
1750  globalCoPointSlavesMapPtr_.reset
1751  (
1752  new mapDistribute
1753  (
1754  move(globalData.map())
1755  )
1756  );
1757 
1758  if (debug)
1759  {
1760  Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
1761  << " finished calculating coupled master to collocated"
1762  << " slave point addressing." << endl;
1763  }
1764 }
1765 
1766 
1767 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1768 
1770 :
1771  processorTopology(mesh.boundaryMesh(), UPstream::worldComm),
1772  mesh_(mesh),
1773  nTotalPoints_(-1),
1774  nTotalFaces_(-1),
1775  nTotalCells_(-1),
1776  processorPatches_(0),
1777  processorPatchIndices_(0),
1778  processorPatchNeighbours_(0),
1779  nGlobalPoints_(-1),
1780  sharedPointLabelsPtr_(nullptr),
1781  sharedPointAddrPtr_(nullptr),
1782  sharedPointGlobalLabelsPtr_(nullptr),
1783  nGlobalEdges_(-1),
1784  sharedEdgeLabelsPtr_(nullptr),
1785  sharedEdgeAddrPtr_(nullptr)
1786 {
1787  updateMesh();
1788 }
1789 
1790 
1791 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1792 
1794 {
1795  clearOut();
1796 }
1797 
1798 
1800 {
1801  // Point
1802  nGlobalPoints_ = -1;
1803  sharedPointLabelsPtr_.clear();
1804  sharedPointAddrPtr_.clear();
1805  sharedPointGlobalLabelsPtr_.clear();
1806 
1807  // Edge
1808  nGlobalEdges_ = -1;
1809  sharedEdgeLabelsPtr_.clear();
1810  sharedEdgeAddrPtr_.clear();
1811 
1812  // Coupled patch
1813  coupledPatchPtr_.clear();
1814  coupledPatchMeshEdgesPtr_.clear();
1815  coupledPatchMeshEdgeMapPtr_.clear();
1816  globalTransformsPtr_.clear();
1817 
1818  // Point
1819  globalPointNumberingPtr_.clear();
1820  globalPointSlavesPtr_.clear();
1821  globalPointTransformedSlavesPtr_.clear();
1822  globalPointSlavesMapPtr_.clear();
1823  // Edge
1824  globalEdgeNumberingPtr_.clear();
1825  globalEdgeSlavesPtr_.clear();
1826  globalEdgeTransformedSlavesPtr_.clear();
1827  globalEdgeOrientationPtr_.clear();
1828  globalEdgeSlavesMapPtr_.clear();
1829 
1830  // Face
1831  globalBoundaryFaceNumberingPtr_.clear();
1832  globalPointBoundaryFacesPtr_.clear();
1833  globalPointTransformedBoundaryFacesPtr_.clear();
1834  globalPointBoundaryFacesMapPtr_.clear();
1835 
1836  // Cell
1837  boundaryCellsPtr_.clear();
1838  globalBoundaryCellNumberingPtr_.clear();
1839  globalPointBoundaryCellsPtr_.clear();
1840  globalPointTransformedBoundaryCellsPtr_.clear();
1841  globalPointBoundaryCellsMapPtr_.clear();
1842 
1843  // Other: collocated points
1844  globalCoPointSlavesPtr_.clear();
1845  globalCoPointSlavesMapPtr_.clear();
1846 }
1847 
1848 
1849 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1850 
1852 {
1853  if (!sharedPointGlobalLabelsPtr_.valid())
1854  {
1855  sharedPointGlobalLabelsPtr_.reset
1856  (
1858  );
1859  labelList& sharedPointGlobalLabels = sharedPointGlobalLabelsPtr_();
1860 
1861  IOobject addrHeader
1862  (
1863  "pointProcAddressing",
1864  mesh_.facesInstance()/mesh_.meshSubDir,
1865  mesh_,
1867  );
1868 
1869  if (addrHeader.typeHeaderOk<labelIOList>(true))
1870  {
1871  // There is a pointProcAddressing file so use it to get labels
1872  // on the original mesh
1873  Pout<< "globalMeshData::sharedPointGlobalLabels : "
1874  << "Reading pointProcAddressing" << endl;
1875 
1876  labelIOList pointProcAddressing(addrHeader);
1877 
1879 
1880  forAll(pointLabels, i)
1881  {
1882  // Get my mesh point
1883  label pointi = pointLabels[i];
1884 
1885  // Map to mesh point of original mesh
1886  sharedPointGlobalLabels[i] = pointProcAddressing[pointi];
1887  }
1888  }
1889  else
1890  {
1891  Pout<< "globalMeshData::sharedPointGlobalLabels :"
1892  << " Setting pointProcAddressing to -1" << endl;
1893 
1894  sharedPointGlobalLabels = -1;
1895  }
1896  }
1897  return sharedPointGlobalLabelsPtr_();
1898 }
1899 
1900 
1902 {
1903  // Get all processors to send their shared points to master.
1904  // (not very efficient)
1905 
1907  const labelList& pointAddr = sharedPointAddr();
1909 
1910  if (Pstream::master())
1911  {
1912  // Master:
1913  // insert my own data first
1914  forAll(pointLabels, i)
1915  {
1916  label sharedPointi = pointAddr[i];
1917 
1918  sharedPoints[sharedPointi] = mesh_.points()[pointLabels[i]];
1919  }
1920 
1921  // Receive data from slaves and insert
1922  for
1923  (
1924  int slave=Pstream::firstSlave();
1925  slave<=Pstream::lastSlave();
1926  slave++
1927  )
1928  {
1929  IPstream fromSlave(Pstream::commsTypes::blocking, slave);
1930 
1931  labelList nbrSharedPointAddr;
1932  pointField nbrSharedPoints;
1933  fromSlave >> nbrSharedPointAddr >> nbrSharedPoints;
1934 
1935  forAll(nbrSharedPointAddr, i)
1936  {
1937  label sharedPointi = nbrSharedPointAddr[i];
1938 
1939  sharedPoints[sharedPointi] = nbrSharedPoints[i];
1940  }
1941  }
1942 
1943  // Send back
1944  for
1945  (
1946  int slave=Pstream::firstSlave();
1947  slave<=Pstream::lastSlave();
1948  slave++
1949  )
1950  {
1951  OPstream toSlave
1952  (
1954  slave,
1955  sharedPoints.size()*sizeof(Zero)
1956  );
1957  toSlave << sharedPoints;
1958  }
1959  }
1960  else
1961  {
1962  // Slave:
1963  // send points
1964  {
1965  OPstream toMaster
1966  (
1969  );
1970  toMaster
1971  << pointAddr
1972  << UIndirectList<point>(mesh_.points(), pointLabels)();
1973  }
1974 
1975  // Receive sharedPoints
1976  {
1977  IPstream fromMaster
1978  (
1981  );
1982  fromMaster >> sharedPoints;
1983  }
1984  }
1985 
1986  return sharedPoints;
1987 }
1988 
1989 
1991 {
1992  // Get coords of my shared points
1994 
1995  // Append from all processors
1997 
1998  // Merge tolerance
1999  scalar tolDim = matchTol_ * mesh_.bounds().mag();
2000 
2001  // And see how many are unique
2002  labelList pMap;
2003  pointField mergedPoints;
2004 
2006  (
2007  sharedPoints, // coordinates to merge
2008  tolDim, // tolerance
2009  false, // verbosity
2010  pMap,
2011  mergedPoints
2012  );
2013 
2014  return mergedPoints;
2015 }
2016 
2017 
2019 {
2020  if (nGlobalPoints_ == -1)
2021  {
2022  calcSharedPoints();
2023  }
2024  return nGlobalPoints_;
2025 }
2026 
2027 
2029 {
2030  if (!sharedPointLabelsPtr_.valid())
2031  {
2032  calcSharedPoints();
2033  }
2034  return sharedPointLabelsPtr_();
2035 }
2036 
2037 
2039 {
2040  if (!sharedPointAddrPtr_.valid())
2041  {
2042  calcSharedPoints();
2043  }
2044  return sharedPointAddrPtr_();
2045 }
2046 
2047 
2049 {
2050  if (nGlobalEdges_ == -1)
2051  {
2052  calcSharedEdges();
2053  }
2054  return nGlobalEdges_;
2055 }
2056 
2057 
2059 {
2060  if (!sharedEdgeLabelsPtr_.valid())
2061  {
2062  calcSharedEdges();
2063  }
2064  return sharedEdgeLabelsPtr_();
2065 }
2066 
2067 
2069 {
2070  if (!sharedEdgeAddrPtr_.valid())
2071  {
2072  calcSharedEdges();
2073  }
2074  return sharedEdgeAddrPtr_();
2075 }
2076 
2077 
2079 {
2080  if (!coupledPatchPtr_.valid())
2081  {
2082  const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
2083 
2084  label nCoupled = 0;
2085 
2086  forAll(bMesh, patchi)
2087  {
2088  const polyPatch& pp = bMesh[patchi];
2089 
2090  if (pp.coupled())
2091  {
2092  nCoupled += pp.size();
2093  }
2094  }
2095  labelList coupledFaces(nCoupled);
2096  nCoupled = 0;
2097 
2098  forAll(bMesh, patchi)
2099  {
2100  const polyPatch& pp = bMesh[patchi];
2101 
2102  if (pp.coupled())
2103  {
2104  label facei = pp.start();
2105 
2106  forAll(pp, i)
2107  {
2108  coupledFaces[nCoupled++] = facei++;
2109  }
2110  }
2111  }
2112 
2113  coupledPatchPtr_.reset
2114  (
2116  (
2118  (
2119  mesh_.faces(),
2120  coupledFaces
2121  ),
2122  mesh_.points()
2123  )
2124  );
2125 
2126  if (debug)
2127  {
2128  Pout<< "globalMeshData::coupledPatch() :"
2129  << " constructed coupled faces patch:"
2130  << " faces:" << coupledPatchPtr_().size()
2131  << " points:" << coupledPatchPtr_().nPoints()
2132  << endl;
2133  }
2134  }
2135  return coupledPatchPtr_();
2136 }
2137 
2138 
2140 {
2141  if (!coupledPatchMeshEdgesPtr_.valid())
2142  {
2143  coupledPatchMeshEdgesPtr_.reset
2144  (
2145  new labelList
2146  (
2147  coupledPatch().meshEdges
2148  (
2149  mesh_.edges(),
2150  mesh_.pointEdges()
2151  )
2152  )
2153  );
2154  }
2155  return coupledPatchMeshEdgesPtr_();
2156 }
2157 
2158 
2160 const
2161 {
2162  if (!coupledPatchMeshEdgeMapPtr_.valid())
2163  {
2164  const labelList& me = coupledPatchMeshEdges();
2165 
2166  coupledPatchMeshEdgeMapPtr_.reset(new Map<label>(2*me.size()));
2167  Map<label>& em = coupledPatchMeshEdgeMapPtr_();
2168 
2169  forAll(me, i)
2170  {
2171  em.insert(me[i], i);
2172  }
2173  }
2174  return coupledPatchMeshEdgeMapPtr_();
2175 }
2176 
2177 
2179 {
2180  if (!globalPointNumberingPtr_.valid())
2181  {
2182  globalPointNumberingPtr_.reset
2183  (
2185  );
2186  }
2187  return globalPointNumberingPtr_();
2188 }
2189 
2190 
2193 {
2194  if (!globalTransformsPtr_.valid())
2195  {
2196  globalTransformsPtr_.reset(new globalIndexAndTransform(mesh_));
2197  }
2198  return globalTransformsPtr_();
2199 }
2200 
2201 
2203 {
2204  if (!globalPointSlavesPtr_.valid())
2205  {
2206  calcGlobalPointSlaves();
2207  }
2208  return globalPointSlavesPtr_();
2209 }
2210 
2211 
2213 const
2214 {
2215  if (!globalPointTransformedSlavesPtr_.valid())
2216  {
2217  calcGlobalPointSlaves();
2218  }
2219  return globalPointTransformedSlavesPtr_();
2220 }
2221 
2222 
2224 {
2225  if (!globalPointSlavesMapPtr_.valid())
2226  {
2227  calcGlobalPointSlaves();
2228  }
2229  return globalPointSlavesMapPtr_();
2230 }
2231 
2232 
2234 {
2235  if (!globalEdgeNumberingPtr_.valid())
2236  {
2237  globalEdgeNumberingPtr_.reset
2238  (
2239  new globalIndex(coupledPatch().nEdges())
2240  );
2241  }
2242  return globalEdgeNumberingPtr_();
2243 }
2244 
2245 
2247 {
2248  if (!globalEdgeSlavesPtr_.valid())
2249  {
2250  calcGlobalEdgeSlaves();
2251  }
2252  return globalEdgeSlavesPtr_();
2253 }
2254 
2255 
2257 const
2258 {
2259  if (!globalEdgeTransformedSlavesPtr_.valid())
2260  {
2261  calcGlobalEdgeSlaves();
2262  }
2263  return globalEdgeTransformedSlavesPtr_();
2264 }
2265 
2266 
2268 {
2269  if (!globalEdgeOrientationPtr_.valid())
2270  {
2271  calcGlobalEdgeOrientation();
2272  }
2273  return globalEdgeOrientationPtr_();
2274 }
2275 
2276 
2278 {
2279  if (!globalEdgeSlavesMapPtr_.valid())
2280  {
2281  calcGlobalEdgeSlaves();
2282  }
2283  return globalEdgeSlavesMapPtr_();
2284 }
2285 
2286 
2288 const
2289 {
2290  if (!globalBoundaryFaceNumberingPtr_.valid())
2291  {
2292  calcGlobalPointBoundaryFaces();
2293  }
2294  return globalBoundaryFaceNumberingPtr_();
2295 }
2296 
2297 
2299 const
2300 {
2301  if (!globalPointBoundaryFacesPtr_.valid())
2302  {
2303  calcGlobalPointBoundaryFaces();
2304  }
2305  return globalPointBoundaryFacesPtr_();
2306 }
2307 
2308 
2309 const Foam::labelListList&
2311 {
2312  if (!globalPointTransformedBoundaryFacesPtr_.valid())
2313  {
2314  calcGlobalPointBoundaryFaces();
2315  }
2316  return globalPointTransformedBoundaryFacesPtr_();
2317 }
2318 
2319 
2321 const
2322 {
2323  if (!globalPointBoundaryFacesMapPtr_.valid())
2324  {
2325  calcGlobalPointBoundaryFaces();
2326  }
2327  return globalPointBoundaryFacesMapPtr_();
2328 }
2329 
2330 
2332 {
2333  if (!boundaryCellsPtr_.valid())
2334  {
2335  calcGlobalPointBoundaryCells();
2336  }
2337  return boundaryCellsPtr_();
2338 }
2339 
2340 
2342 const
2343 {
2344  if (!globalBoundaryCellNumberingPtr_.valid())
2345  {
2346  calcGlobalPointBoundaryCells();
2347  }
2348  return globalBoundaryCellNumberingPtr_();
2349 }
2350 
2351 
2353 const
2354 {
2355  if (!globalPointBoundaryCellsPtr_.valid())
2356  {
2357  calcGlobalPointBoundaryCells();
2358  }
2359  return globalPointBoundaryCellsPtr_();
2360 }
2361 
2362 
2363 const Foam::labelListList&
2365 {
2366  if (!globalPointTransformedBoundaryCellsPtr_.valid())
2367  {
2368  calcGlobalPointBoundaryCells();
2369  }
2370  return globalPointTransformedBoundaryCellsPtr_();
2371 }
2372 
2373 
2375 const
2376 {
2377  if (!globalPointBoundaryCellsMapPtr_.valid())
2378  {
2379  calcGlobalPointBoundaryCells();
2380  }
2381  return globalPointBoundaryCellsMapPtr_();
2382 }
2383 
2384 
2386 {
2387  if (!globalCoPointSlavesPtr_.valid())
2388  {
2389  calcGlobalCoPointSlaves();
2390  }
2391  return globalCoPointSlavesPtr_();
2392 }
2393 
2394 
2396 {
2397  if (!globalCoPointSlavesMapPtr_.valid())
2398  {
2399  calcGlobalCoPointSlaves();
2400  }
2401  return globalCoPointSlavesMapPtr_();
2402 }
2403 
2404 
2407  labelList& pointToGlobal,
2408  labelList& uniquePoints
2409 ) const
2410 {
2411  const indirectPrimitivePatch& cpp = coupledPatch();
2412  const globalIndex& globalCoupledPoints = globalPointNumbering();
2413  // Use collocated only
2414  const labelListList& pointSlaves = globalCoPointSlaves();
2415  const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
2416 
2417 
2418  // Points are either
2419  // - master with slaves
2420  // - slave with a master
2421  // - other (since e.g. non-collocated cyclics not connected)
2422 
2423  labelList masterGlobalPoint(cpp.nPoints(), -1);
2424  forAll(masterGlobalPoint, pointi)
2425  {
2426  const labelList& slavePoints = pointSlaves[pointi];
2427  if (slavePoints.size() > 0)
2428  {
2429  masterGlobalPoint[pointi] = globalCoupledPoints.toGlobal(pointi);
2430  }
2431  }
2432 
2433  // Sync by taking max
2434  syncData
2435  (
2436  masterGlobalPoint,
2437  pointSlaves,
2438  labelListList(0), // no transforms
2439  pointSlavesMap,
2440  maxEqOp<label>()
2441  );
2442 
2443 
2444  // 1. Count number of masters on my processor.
2445  label nMaster = 0;
2446  PackedBoolList isMaster(mesh_.nPoints(), 1);
2447  forAll(pointSlaves, pointi)
2448  {
2449  if (masterGlobalPoint[pointi] == -1)
2450  {
2451  // unconnected point (e.g. from separated cyclic)
2452  nMaster++;
2453  }
2454  else if
2455  (
2456  masterGlobalPoint[pointi]
2457  == globalCoupledPoints.toGlobal(pointi)
2458  )
2459  {
2460  // connected master
2461  nMaster++;
2462  }
2463  else
2464  {
2465  // connected slave point
2466  isMaster[cpp.meshPoints()[pointi]] = 0;
2467  }
2468  }
2469 
2470  label myUniquePoints = mesh_.nPoints() - cpp.nPoints() + nMaster;
2471 
2472  // Pout<< "Points :" << nl
2473  // << " mesh : " << mesh_.nPoints() << nl
2474  // << " of which coupled : " << cpp.nPoints() << nl
2475  // << " of which master : " << nMaster << nl
2476  // << endl;
2477 
2478 
2479  // 2. Create global indexing for unique points.
2480  autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
2481 
2482 
2483  // 3. Assign global point numbers. Keep slaves unset.
2484  pointToGlobal.setSize(mesh_.nPoints());
2485  pointToGlobal = -1;
2486  uniquePoints.setSize(myUniquePoints);
2487  nMaster = 0;
2488 
2489  forAll(isMaster, meshPointi)
2490  {
2491  if (isMaster[meshPointi])
2492  {
2493  pointToGlobal[meshPointi] = globalPointsPtr().toGlobal(nMaster);
2494  uniquePoints[nMaster] = meshPointi;
2495  nMaster++;
2496  }
2497  }
2498 
2499 
2500  // 4. Push global index for coupled points to slaves.
2501  {
2502  labelList masterToGlobal(pointSlavesMap.constructSize(), -1);
2503 
2504  forAll(pointSlaves, pointi)
2505  {
2506  const labelList& slaves = pointSlaves[pointi];
2507 
2508  if (slaves.size() > 0)
2509  {
2510  // Duplicate master globalpoint into slave slots
2511  label meshPointi = cpp.meshPoints()[pointi];
2512  masterToGlobal[pointi] = pointToGlobal[meshPointi];
2513  forAll(slaves, i)
2514  {
2515  masterToGlobal[slaves[i]] = masterToGlobal[pointi];
2516  }
2517  }
2518  }
2519 
2520  // Send back
2521  pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
2522 
2523  // On slave copy master index into overal map.
2524  forAll(pointSlaves, pointi)
2525  {
2526  label meshPointi = cpp.meshPoints()[pointi];
2527 
2528  if (!isMaster[meshPointi])
2529  {
2530  pointToGlobal[meshPointi] = masterToGlobal[pointi];
2531  }
2532  }
2533  }
2534 
2535  return globalPointsPtr;
2536 }
2537 
2538 
2541  const labelList& meshPoints,
2542  const Map<label>& meshPointMap,
2543  labelList& pointToGlobal,
2544  labelList& uniqueMeshPoints
2545 ) const
2546 {
2547  const indirectPrimitivePatch& cpp = coupledPatch();
2548  const labelListList& pointSlaves = globalCoPointSlaves();
2549  const mapDistribute& pointSlavesMap = globalCoPointSlavesMap();
2550 
2551 
2552  // The patch points come in two variants:
2553  // - not on a coupled patch so guaranteed unique
2554  // - on a coupled patch
2555  // If the point is on a coupled patch the problem is that the
2556  // master-slave structure (globalPointSlaves etc.) assigns one of the
2557  // coupled points to be the master but this master point is not
2558  // necessarily on the patch itself! (it might just be connected to the
2559  // patch point via coupled patches).
2560 
2561 
2562  // Determine mapping:
2563  // - from patch point to coupled point (or -1)
2564  // - from coupled point to global patch point
2565  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2566 
2567  globalIndex globalPPoints(meshPoints.size());
2568 
2569  labelList patchToCoupled(meshPoints.size(), -1);
2570  label nCoupled = 0;
2571  labelList coupledToGlobalPatch(pointSlavesMap.constructSize(), -1);
2572 
2573  // Note: loop over patch since usually smaller
2574  forAll(meshPoints, patchPointi)
2575  {
2576  label meshPointi = meshPoints[patchPointi];
2577 
2578  Map<label>::const_iterator iter = cpp.meshPointMap().find(meshPointi);
2579 
2580  if (iter != cpp.meshPointMap().end())
2581  {
2582  patchToCoupled[patchPointi] = iter();
2583  coupledToGlobalPatch[iter()] = globalPPoints.toGlobal(patchPointi);
2584  nCoupled++;
2585  }
2586  }
2587 
2588  // Pout<< "Patch:" << nl
2589  // << " points:" << meshPoints.size() << nl
2590  // << " of which on coupled patch:" << nCoupled << endl;
2591 
2592 
2593  // Determine master of connected points
2594  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2595  // Problem is that the coupled master might not be on the patch. So
2596  // work out the best patch-point master from all connected points.
2597  // - if the coupled master is on the patch it becomes the patch-point master
2598  // - else the slave with the lowest numbered patch point label
2599 
2600  // Get all data on master
2601  pointSlavesMap.distribute(coupledToGlobalPatch);
2602  forAll(pointSlaves, coupledPointi)
2603  {
2604  const labelList& slaves = pointSlaves[coupledPointi];
2605 
2606  if (slaves.size() > 0)
2607  {
2608  // I am master. What is the best candidate for patch-point master
2609  label masterI = labelMax;
2610  if (coupledToGlobalPatch[coupledPointi] != -1)
2611  {
2612  // I am master and on the coupled patch. Use me.
2613  masterI = coupledToGlobalPatch[coupledPointi];
2614  }
2615  else
2616  {
2617  // Get min of slaves as master.
2618  forAll(slaves, i)
2619  {
2620  label slavePp = coupledToGlobalPatch[slaves[i]];
2621  if (slavePp != -1 && slavePp < masterI)
2622  {
2623  masterI = slavePp;
2624  }
2625  }
2626  }
2627 
2628  if (masterI != labelMax)
2629  {
2630  // Push back
2631  coupledToGlobalPatch[coupledPointi] = masterI;
2632  forAll(slaves, i)
2633  {
2634  coupledToGlobalPatch[slaves[i]] = masterI;
2635  }
2636  }
2637  }
2638  }
2639  pointSlavesMap.reverseDistribute(cpp.nPoints(), coupledToGlobalPatch);
2640 
2641 
2642  // Generate compact numbering for master points
2643  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2644  // Now coupledToGlobalPatch is the globalIndex of the master point.
2645  // Now every processor can check whether they hold it and generate a
2646  // compact numbering.
2647 
2648  label nMasters = 0;
2649  forAll(meshPoints, patchPointi)
2650  {
2651  if (patchToCoupled[patchPointi] == -1)
2652  {
2653  nMasters++;
2654  }
2655  else
2656  {
2657  label coupledPointi = patchToCoupled[patchPointi];
2658  if
2659  (
2660  globalPPoints.toGlobal(patchPointi)
2661  == coupledToGlobalPatch[coupledPointi]
2662  )
2663  {
2664  // I am the master
2665  nMasters++;
2666  }
2667  }
2668  }
2669 
2670  autoPtr<globalIndex> globalPointsPtr(new globalIndex(nMasters));
2671 
2672  // Pout<< "Patch:" << nl
2673  // << " points:" << meshPoints.size() << nl
2674  // << " of which on coupled patch:" << nCoupled << nl
2675  // << " of which master:" << nMasters << endl;
2676 
2677 
2678 
2679  // Push back compact numbering for master point onto slaves
2680  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2681 
2682  pointToGlobal.setSize(meshPoints.size());
2683  pointToGlobal = -1;
2684  uniqueMeshPoints.setSize(nMasters);
2685 
2686  // Sync master in global point numbering so all know the master point.
2687  // Initialise globalMaster to be -1 except at a globalMaster.
2688  labelList globalMaster(cpp.nPoints(), -1);
2689 
2690  nMasters = 0;
2691  forAll(meshPoints, patchPointi)
2692  {
2693  if (patchToCoupled[patchPointi] == -1)
2694  {
2695  uniqueMeshPoints[nMasters++] = meshPoints[patchPointi];
2696  }
2697  else
2698  {
2699  label coupledPointi = patchToCoupled[patchPointi];
2700  if
2701  (
2702  globalPPoints.toGlobal(patchPointi)
2703  == coupledToGlobalPatch[coupledPointi]
2704  )
2705  {
2706  globalMaster[coupledPointi] =
2707  globalPointsPtr().toGlobal(nMasters);
2708  uniqueMeshPoints[nMasters++] = meshPoints[patchPointi];
2709  }
2710  }
2711  }
2712 
2713 
2714  // Sync by taking max
2715  syncData
2716  (
2717  globalMaster,
2718  pointSlaves,
2719  labelListList(0), // no transforms
2720  pointSlavesMap,
2721  maxEqOp<label>()
2722  );
2723 
2724 
2725  // Now everyone has the master point in globalPointsPtr numbering. Fill
2726  // in the pointToGlobal map.
2727  nMasters = 0;
2728  forAll(meshPoints, patchPointi)
2729  {
2730  if (patchToCoupled[patchPointi] == -1)
2731  {
2732  pointToGlobal[patchPointi] = globalPointsPtr().toGlobal(nMasters++);
2733  }
2734  else
2735  {
2736  label coupledPointi = patchToCoupled[patchPointi];
2737  pointToGlobal[patchPointi] = globalMaster[coupledPointi];
2738 
2739  if
2740  (
2741  globalPPoints.toGlobal(patchPointi)
2742  == coupledToGlobalPatch[coupledPointi]
2743  )
2744  {
2745  nMasters++;
2746  }
2747  }
2748  }
2749 
2750  return globalPointsPtr;
2751 }
2752 
2753 
2755 {
2756  // Topology does not change and we don't store any geometry so nothing
2757  // needs to be done.
2758  // Only global transformations might change but this is not really
2759  // supported.
2760 }
2761 
2762 
2764 {
2765  // Clear out old data
2766  clearOut();
2767 
2768  // Do processor patch addressing
2769  initProcAddr();
2770 
2771  scalar tolDim = matchTol_ * mesh_.bounds().mag();
2772 
2773  if (debug)
2774  {
2775  Pout<< "globalMeshData : merge dist:" << tolDim << endl;
2776  }
2777 
2778  // *** Temporary hack to avoid problems with overlapping communication
2779  // *** between these reductions and the calculation of deltaCoeffs
2781  (
2784  true
2785  );
2786 
2787  // Total number of faces.
2788  nTotalFaces_ = returnReduce
2789  (
2790  mesh_.nFaces(),
2791  sumOp<label>(),
2792  Pstream::msgType(),
2793  comm
2794  );
2795 
2796  if (debug)
2797  {
2798  Pout<< "globalMeshData : nTotalFaces_:" << nTotalFaces_ << endl;
2799  }
2800 
2801  nTotalCells_ = returnReduce
2802  (
2803  mesh_.nCells(),
2804  sumOp<label>(),
2805  Pstream::msgType(),
2806  comm
2807  );
2808 
2809  if (debug)
2810  {
2811  Pout<< "globalMeshData : nTotalCells_:" << nTotalCells_ << endl;
2812  }
2813 
2814  nTotalPoints_ = returnReduce
2815  (
2816  mesh_.nPoints(),
2817  sumOp<label>(),
2818  Pstream::msgType(),
2819  comm
2820  );
2821 
2823 
2824  if (debug)
2825  {
2826  Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;
2827  }
2828 }
2829 
2830 
2831 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
const labelList & sharedPointLabels() const
Return indices of local points that are globally shared.
const Map< label > & coupledPatchMeshEdgeMap() const
Return map from mesh edges to coupledPatch edges.
transformer transforms
Less function class used in sorting encoded transforms and indices.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
label nPoints() const
Return number of points supporting patch faces.
const labelListList & transformedPointPoints() const
Transformed points per point (in mapDistribute indices)
Definition: globalPoints.H:332
const labelListList & globalPointSlaves() const
const mapDistribute & globalCoPointSlavesMap() const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void combineReduce(const List< UPstream::commsStruct > &comms, T &Value, const CombineOp &cop, const int tag, const label comm)
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
static int masterNo()
Process index of the master.
Definition: UPstream.H:417
label index(const labelPair &globalIAndTransform) const
Index carried by the object.
const labelListList & globalPointTransformedSlaves() const
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
pointField geometricSharedPoints() const
Like sharedPoints but keeps cyclic points separate.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:808
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
const globalIndex & globalBoundaryCellNumbering() const
Numbering of boundary cells is according to boundaryCells()
labelList pointLabels(nPoints, -1)
const labelList & sharedPointGlobalLabels() const
Return shared point global labels. Tries to read.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
const labelListList & globalPointTransformedBoundaryFaces() const
static int firstSlave()
Process index of first slave.
Definition: UPstream.H:446
const labelListList & pointEdges() const
const mapDistribute & globalEdgeSlavesMap() const
label nFaces() const
static const Foam::scalar matchTol_
Geometric tolerance (fraction of bounding box)
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:312
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int compare(const Pair< Type > &a, const Pair< Type > &b)
Compare Pairs.
Definition: Pair.H:153
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
void clearOut()
Remove all demand driven data.
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const mapDistribute & globalPointSlavesMap() const
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
const labelListList & globalPointBoundaryCells() const
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
label processor(const labelPair &globalIAndTransform) const
Which processor does this come from?
static int & msgType()
Message tag of standard messages.
Definition: UPstream.H:476
autoPtr< globalIndex > mergePoints(labelList &pointToGlobal, labelList &uniquePoints) const
Helper for merging (collocated!) mesh point data.
const labelListList & globalEdgeSlaves() const
static label worldComm
Default communicator (all processors)
Definition: UPstream.H:278
const labelList & coupledPatchMeshEdges() const
Return map from coupledPatch edges to mesh edges.
const labelListList & globalPointBoundaryFaces() const
const mapDistribute & map() const
Corresponding map.
Definition: globalPoints.H:344
labelPair encode(const label index, const label transformIndex) const
Encode index and bare index as components on own processor.
Combination-Reduction operation for a parallel run. The information from all nodes is collected on th...
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Determines processor-processor connection. After instantiation contains on all processors the process...
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
Input inter-processor communications stream.
Definition: IPstream.H:50
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
const labelList & boundaryCells() const
From boundary cell to mesh cell.
label subtractTransformIndex(const label transformIndex0, const label transformIndex1) const
Subtract two transformIndices.
label nGlobalEdges() const
Return number of globally shared edges. Demand-driven.
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
scalar y
A list of faces which address into the list of points.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
A List obtained as a section of another List.
Definition: SubList.H:53
const globalIndex & globalEdgeNumbering() const
label nGlobalPoints() const
Return number of globally shared points.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
const labelList & sharedEdgeAddr() const
Return addressing into the complete globally shared edge.
dynamicFvMesh & mesh
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:319
const globalIndexAndTransform & globalTransforms() const
Global transforms numbering.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
const PackedBoolList & globalEdgeOrientation() const
Is my edge same orientation as master edge.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
label nPoints
void clear()
Clear all entries from table.
Definition: HashTable.C:468
void sort(UList< T > &)
Definition: UList.C:115
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:90
const globalIndex & globalPointNumbering() const
Numbering of coupled points is according to coupledPatch.
label size() const
Global sum of localSizes.
Definition: globalIndexI.H:66
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
label nullTransformIndex() const
Return the transformIndex (index in transformPermutations)
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
label toLocal(const label i) const
From global to local on current processor.
Definition: globalIndexI.H:117
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
label whichProcID(const label i) const
Which processor does global come from? Binary search.
Definition: globalIndexI.H:123
List< labelPair > labelPairList
List of labelPairs.
Definition: labelPair.H:49
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
static const char nl
Definition: Ostream.H:260
Map from edge (expressed as its endpoints) to value.
Definition: EdgeMap.H:47
const labelListList & pointPoints() const
Non-transformed connected points per point (in mapDistribute.
Definition: globalPoints.H:319
~globalMeshData()
Destructor.
defineTypeNameAndDebug(combustionModel, 0)
Merge points. See below.
labelList f(nPoints)
Output inter-processor communications stream.
Definition: OPstream.H:50
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void movePoints(const pointField &newPoints)
Update for moving points.
label transformIndex(const labelPair &globalIAndTransform) const
Transform carried by the object.
void reverseDistribute(const label constructSize, List< T > &, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Reverse distribute data using default commsType.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
static void syncData(List< Type > &pointData, const labelListList &slaves, const labelListList &transformedSlaves, const mapDistribute &slavesMap, const globalIndexAndTransform &, const CombineOp &cop, const TransformOp &top)
Helper: synchronise data with transforms.
label constructSize() const
Constructed data size.
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:100
void setSize(const label)
Reset size of List.
Definition: List.C:281
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:440
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
const labelListList & globalPointTransformedBoundaryCells() const
A bit-packed bool list.
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
label patchi
Class containing processor-to-processor mapping information.
const dimensionedScalar me
Electron mass.
globalMeshData(const polyMesh &mesh)
Construct from mesh, derive rest (does parallel communication!)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const polyBoundaryMesh & bMesh
A List with indirect addressing.
Definition: fvMatrix.H:106
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
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.
void updateMesh()
Change global mesh data given a topological change. Does a.
label n
label nPoints() const
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static label allocateCommunicator(const label parent, const labelList &subRanks, const bool doPstream=true)
Allocate a new communicator.
Definition: UPstream.C:248
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
const mapDistribute & globalPointBoundaryFacesMap() const
const labelList & sharedPointAddr() const
Return addressing into the complete globally shared points.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
pointField sharedPoints() const
Collect coordinates of shared points on all processors.
const mapDistribute & globalPointBoundaryCellsMap() const
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
A List with indirect addressing.
Definition: IndirectList.H:101
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
Inter-processor communications stream.
Definition: UPstream.H:58
static void freeCommunicator(const label communicator, const bool doPstream=true)
Free a previously allocated communicator.
Definition: UPstream.C:314
const labelList & sharedEdgeLabels() const
Return indices of local edges that are globally shared.
Namespace for OpenFOAM.
Determination and storage of the possible independent transforms introduced by coupledPolyPatches, as well as all of the possible permutations of these transforms generated by the presence of multiple coupledPolyPatches, i.e. more than one cyclic boundary. Note that any given point can be on maximum 3 transforms only (and these transforms have to be perpendicular)
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477
const labelListList & globalCoPointSlaves() const
static int lastSlave(const label communicator=0)
Process index of last slave.
Definition: UPstream.H:452
const labelListList & globalEdgeTransformedSlaves() const
void operator()(labelPair &x, const labelPair &y) const
label localSize() const
My local size.
Definition: globalIndexI.H:60
const globalIndex & globalBoundaryFaceNumbering() const
Numbering of boundary faces is face-mesh.nInternalFaces()