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