PointEdgeWave.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 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 "PointEdgeWave.H"
27 #include "polyMesh.H"
28 #include "processorPolyPatch.H"
29 #include "cyclicPolyPatch.H"
30 #include "OPstream.H"
31 #include "IPstream.H"
33 #include "debug.H"
34 #include "typeInfo.H"
35 #include "globalMeshData.H"
36 #include "pointFields.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 template<class Type, class TrackingData>
42 
43 template<class Type, class TrackingData>
45 
46 namespace Foam
47 {
48  //- Reduction class. If x and y are not equal assign value.
49  template<class Type, class TrackingData>
51  {
52  TrackingData& td_;
53 
54  public:
55  combineEqOp(TrackingData& td)
56  :
57  td_(td)
58  {}
59 
60  void operator()(Type& x, const Type& y) const
61  {
62  if (!x.valid(td_) && y.valid(td_))
63  {
64  x = y;
65  }
66  }
67  };
68 }
69 
70 
71 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
72 
73 // Handle leaving domain. Implementation referred to Type
74 template<class Type, class TrackingData>
76 (
77  const polyPatch& patch,
78  const labelList& patchPointLabels,
79  List<Type>& pointInfo
80 ) const
81 {
82  const labelList& meshPoints = patch.meshPoints();
83 
84  forAll(patchPointLabels, i)
85  {
86  label patchPointI = patchPointLabels[i];
87 
88  const point& pt = patch.points()[meshPoints[patchPointI]];
89 
90  pointInfo[i].leaveDomain(patch, patchPointI, pt, td_);
91  }
92 }
93 
94 
95 // Handle entering domain. Implementation referred to Type
96 template<class Type, class TrackingData>
98 (
99  const polyPatch& patch,
100  const labelList& patchPointLabels,
101  List<Type>& pointInfo
102 ) const
103 {
104  const labelList& meshPoints = patch.meshPoints();
105 
106  forAll(patchPointLabels, i)
107  {
108  label patchPointI = patchPointLabels[i];
109 
110  const point& pt = patch.points()[meshPoints[patchPointI]];
111 
112  pointInfo[i].enterDomain(patch, patchPointI, pt, td_);
113  }
114 }
115 
116 
117 // Transform. Implementation referred to Type
118 template<class Type, class TrackingData>
120 (
121  const polyPatch& patch,
122  const tensorField& rotTensor,
123  List<Type>& pointInfo
124 ) const
125 {
126  if (rotTensor.size() == 1)
127  {
128  const tensor& T = rotTensor[0];
129 
130  forAll(pointInfo, i)
131  {
132  pointInfo[i].transform(T, td_);
133  }
134  }
135  else
136  {
138  (
139  "PointEdgeWave<Type, TrackingData>::transform"
140  "(const tensorField&, List<Type>&)"
141  ) << "Non-uniform transformation on patch " << patch.name()
142  << " of type " << patch.type()
143  << " not supported for point fields"
144  << abort(FatalError);
145 
146  forAll(pointInfo, i)
147  {
148  pointInfo[i].transform(rotTensor[i], td_);
149  }
150  }
151 }
152 
153 
154 // Update info for pointI, at position pt, with information from
155 // neighbouring edge.
156 // Updates:
157 // - changedPoint_, changedPoints_, nChangedPoints_,
158 // - statistics: nEvals_, nUnvisitedPoints_
159 template<class Type, class TrackingData>
161 (
162  const label pointI,
163  const label neighbourEdgeI,
164  const Type& neighbourInfo,
165  Type& pointInfo
166 )
167 {
168  nEvals_++;
169 
170  bool wasValid = pointInfo.valid(td_);
171 
172  bool propagate =
173  pointInfo.updatePoint
174  (
175  mesh_,
176  pointI,
177  neighbourEdgeI,
178  neighbourInfo,
179  propagationTol_,
180  td_
181  );
182 
183  if (propagate)
184  {
185  if (!changedPoint_[pointI])
186  {
187  changedPoint_[pointI] = true;
188  changedPoints_[nChangedPoints_++] = pointI;
189  }
190  }
191 
192  if (!wasValid && pointInfo.valid(td_))
193  {
194  --nUnvisitedPoints_;
195  }
196 
197  return propagate;
198 }
199 
200 
201 // Update info for pointI, at position pt, with information from
202 // same point.
203 // Updates:
204 // - changedPoint_, changedPoints_, nChangedPoints_,
205 // - statistics: nEvals_, nUnvisitedPoints_
206 template<class Type, class TrackingData>
208 (
209  const label pointI,
210  const Type& neighbourInfo,
211  Type& pointInfo
212 )
213 {
214  nEvals_++;
215 
216  bool wasValid = pointInfo.valid(td_);
217 
218  bool propagate =
219  pointInfo.updatePoint
220  (
221  mesh_,
222  pointI,
223  neighbourInfo,
224  propagationTol_,
225  td_
226  );
227 
228  if (propagate)
229  {
230  if (!changedPoint_[pointI])
231  {
232  changedPoint_[pointI] = true;
233  changedPoints_[nChangedPoints_++] = pointI;
234  }
235  }
236 
237  if (!wasValid && pointInfo.valid(td_))
238  {
239  --nUnvisitedPoints_;
240  }
241 
242  return propagate;
243 }
244 
245 
246 // Update info for edgeI, at position pt, with information from
247 // neighbouring point.
248 // Updates:
249 // - changedEdge_, changedEdges_, nChangedEdges_,
250 // - statistics: nEvals_, nUnvisitedEdge_
251 template<class Type, class TrackingData>
253 (
254  const label edgeI,
255  const label neighbourPointI,
256  const Type& neighbourInfo,
257  Type& edgeInfo
258 )
259 {
260  nEvals_++;
261 
262  bool wasValid = edgeInfo.valid(td_);
263 
264  bool propagate =
265  edgeInfo.updateEdge
266  (
267  mesh_,
268  edgeI,
269  neighbourPointI,
270  neighbourInfo,
271  propagationTol_,
272  td_
273  );
274 
275  if (propagate)
276  {
277  if (!changedEdge_[edgeI])
278  {
279  changedEdge_[edgeI] = true;
280  changedEdges_[nChangedEdges_++] = edgeI;
281  }
282  }
283 
284  if (!wasValid && edgeInfo.valid(td_))
285  {
286  --nUnvisitedEdges_;
287  }
288 
289  return propagate;
290 }
291 
292 
293 // Check if patches of given type name are present
294 template<class Type, class TrackingData>
295 template<class PatchType>
297 {
298  label nPatches = 0;
299 
300  forAll(mesh_.boundaryMesh(), patchI)
301  {
302  if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
303  {
304  nPatches++;
305  }
306  }
307  return nPatches;
308 }
309 
310 
311 // Transfer all the information to/from neighbouring processors
312 template<class Type, class TrackingData>
314 {
315  // 1. Send all point info on processor patches.
316 
318 
319  DynamicList<Type> patchInfo;
320  DynamicList<label> thisPoints;
321  DynamicList<label> nbrPoints;
322 
323  forAll(mesh_.globalData().processorPatches(), i)
324  {
325  label patchI = mesh_.globalData().processorPatches()[i];
326  const processorPolyPatch& procPatch =
327  refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
328 
329  patchInfo.clear();
330  patchInfo.reserve(procPatch.nPoints());
331  thisPoints.clear();
332  thisPoints.reserve(procPatch.nPoints());
333  nbrPoints.clear();
334  nbrPoints.reserve(procPatch.nPoints());
335 
336  // Get all changed points in reverse order
337  const labelList& neighbPoints = procPatch.neighbPoints();
338  forAll(neighbPoints, thisPointI)
339  {
340  label meshPointI = procPatch.meshPoints()[thisPointI];
341  if (changedPoint_[meshPointI])
342  {
343  patchInfo.append(allPointInfo_[meshPointI]);
344  thisPoints.append(thisPointI);
345  nbrPoints.append(neighbPoints[thisPointI]);
346  }
347  }
348 
349  // Adapt for leaving domain
350  leaveDomain(procPatch, thisPoints, patchInfo);
351 
352  //if (debug)
353  //{
354  // Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
355  // << " communicating with " << procPatch.neighbProcNo()
356  // << " Sending:" << patchInfo.size() << endl;
357  //}
358 
359  UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
360  toNeighbour << nbrPoints << patchInfo;
361  }
362 
363 
364  pBufs.finishedSends();
365 
366  //
367  // 2. Receive all point info on processor patches.
368  //
369 
370  forAll(mesh_.globalData().processorPatches(), i)
371  {
372  label patchI = mesh_.globalData().processorPatches()[i];
373  const processorPolyPatch& procPatch =
374  refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
375 
376  List<Type> patchInfo;
377  labelList patchPoints;
378 
379  {
380  UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
381  fromNeighbour >> patchPoints >> patchInfo;
382  }
383 
384  //if (debug)
385  //{
386  // Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
387  // << " communicating with " << procPatch.neighbProcNo()
388  // << " Received:" << patchInfo.size() << endl;
389  //}
390 
391  // Apply transform to received data for non-parallel planes
392  if (!procPatch.parallel())
393  {
394  transform(procPatch, procPatch.forwardT(), patchInfo);
395  }
396 
397  // Adapt for entering domain
398  enterDomain(procPatch, patchPoints, patchInfo);
399 
400  // Merge received info
401  const labelList& meshPoints = procPatch.meshPoints();
402  forAll(patchInfo, i)
403  {
404  label meshPointI = meshPoints[patchPoints[i]];
405 
406  if (!allPointInfo_[meshPointI].equal(patchInfo[i], td_))
407  {
408  updatePoint
409  (
410  meshPointI,
411  patchInfo[i],
412  allPointInfo_[meshPointI]
413  );
414  }
415  }
416  }
417 
418  // Collocated points should be handled by face based transfer
419  // (since that is how connectivity is worked out)
420  // They are also explicitly equalised in handleCollocatedPoints to
421  // guarantee identical values.
422 }
423 
424 
425 template<class Type, class TrackingData>
427 {
428  // 1. Send all point info on cyclic patches.
429 
430  DynamicList<Type> nbrInfo;
431  DynamicList<label> nbrPoints;
432  DynamicList<label> thisPoints;
433 
434  forAll(mesh_.boundaryMesh(), patchI)
435  {
436  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
437 
438  if (isA<cyclicPolyPatch>(patch))
439  {
440  const cyclicPolyPatch& cycPatch =
441  refCast<const cyclicPolyPatch>(patch);
442 
443  nbrInfo.clear();
444  nbrInfo.reserve(cycPatch.nPoints());
445  nbrPoints.clear();
446  nbrPoints.reserve(cycPatch.nPoints());
447  thisPoints.clear();
448  thisPoints.reserve(cycPatch.nPoints());
449 
450  // Collect nbrPatch points that have changed
451  {
452  const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
453  const edgeList& pairs = cycPatch.coupledPoints();
454  const labelList& meshPoints = nbrPatch.meshPoints();
455 
456  forAll(pairs, pairI)
457  {
458  label thisPointI = pairs[pairI][0];
459  label nbrPointI = pairs[pairI][1];
460  label meshPointI = meshPoints[nbrPointI];
461 
462  if (changedPoint_[meshPointI])
463  {
464  nbrInfo.append(allPointInfo_[meshPointI]);
465  nbrPoints.append(nbrPointI);
466  thisPoints.append(thisPointI);
467  }
468  }
469 
470  // nbr : adapt for leaving domain
471  leaveDomain(nbrPatch, nbrPoints, nbrInfo);
472  }
473 
474 
475  // Apply rotation for non-parallel planes
476 
477  if (!cycPatch.parallel())
478  {
479  // received data from half1
480  transform(cycPatch, cycPatch.forwardT(), nbrInfo);
481  }
482 
483  //if (debug)
484  //{
485  // Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
486  // << " Changed : " << nbrInfo.size()
487  // << endl;
488  //}
489 
490  // Adapt for entering domain
491  enterDomain(cycPatch, thisPoints, nbrInfo);
492 
493  // Merge received info
494  const labelList& meshPoints = cycPatch.meshPoints();
495  forAll(nbrInfo, i)
496  {
497  label meshPointI = meshPoints[thisPoints[i]];
498 
499  if (!allPointInfo_[meshPointI].equal(nbrInfo[i], td_))
500  {
501  updatePoint
502  (
503  meshPointI,
504  nbrInfo[i],
505  allPointInfo_[meshPointI]
506  );
507  }
508  }
509  }
510  }
511 }
512 
513 
514 // Guarantee collocated points have same information.
515 // Return number of points changed.
516 template<class Type, class TrackingData>
518 {
519  // Transfer onto coupled patch
520  const globalMeshData& gmd = mesh_.globalData();
521  const indirectPrimitivePatch& cpp = gmd.coupledPatch();
522  const labelList& meshPoints = cpp.meshPoints();
523 
524  const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
525  const labelListList& slaves = gmd.globalPointSlaves();
526 
527  List<Type> elems(slavesMap.constructSize());
528  forAll(meshPoints, pointI)
529  {
530  elems[pointI] = allPointInfo_[meshPoints[pointI]];
531  }
532 
533  // Pull slave data onto master (which might or might not have any
534  // initialised points). No need to update transformed slots.
535  slavesMap.distribute(elems, false);
536 
537  // Combine master data with slave data
539 
540  forAll(slaves, pointI)
541  {
542  Type& elem = elems[pointI];
543 
544  const labelList& slavePoints = slaves[pointI];
545 
546  // Combine master with untransformed slave data
547  forAll(slavePoints, j)
548  {
549  cop(elem, elems[slavePoints[j]]);
550  }
551 
552  // Copy result back to slave slots
553  forAll(slavePoints, j)
554  {
555  elems[slavePoints[j]] = elem;
556  }
557  }
558 
559  // Push slave-slot data back to slaves
560  slavesMap.reverseDistribute(elems.size(), elems, false);
561 
562  // Extract back onto mesh
563  forAll(meshPoints, pointI)
564  {
565  if (elems[pointI].valid(td_))
566  {
567  label meshPointI = meshPoints[pointI];
568 
569  Type& elem = allPointInfo_[meshPointI];
570 
571  bool wasValid = elem.valid(td_);
572 
573  // Like updatePoint but bypass Type::updatePoint with its tolerance
574  // checking
575  //if (!elem.valid(td_) || !elem.equal(elems[pointI], td_))
576  if (!elem.equal(elems[pointI], td_))
577  {
578  nEvals_++;
579  elem = elems[pointI];
580 
581  // See if element now valid
582  if (!wasValid && elem.valid(td_))
583  {
584  --nUnvisitedPoints_;
585  }
586 
587  // Update database of changed points
588  if (!changedPoint_[meshPointI])
589  {
590  changedPoint_[meshPointI] = true;
591  changedPoints_[nChangedPoints_++] = meshPointI;
592  }
593  }
594  }
595  }
596 
597  // Sum nChangedPoints over all procs
598  label totNChanged = nChangedPoints_;
599 
600  reduce(totNChanged, sumOp<label>());
601 
602  return totNChanged;
603 }
604 
605 
606 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
607 
608 // Iterate, propagating changedPointsInfo across mesh, until no change (or
609 // maxIter reached). Initial point values specified.
610 template<class Type, class TrackingData>
612 (
613  const polyMesh& mesh,
614  const labelList& changedPoints,
615  const List<Type>& changedPointsInfo,
616 
617  UList<Type>& allPointInfo,
618  UList<Type>& allEdgeInfo,
619 
620  const label maxIter,
621  TrackingData& td
622 )
623 :
624  mesh_(mesh),
625  allPointInfo_(allPointInfo),
626  allEdgeInfo_(allEdgeInfo),
627  td_(td),
628  changedPoint_(mesh_.nPoints(), false),
629  changedPoints_(mesh_.nPoints()),
630  nChangedPoints_(0),
631  changedEdge_(mesh_.nEdges(), false),
632  changedEdges_(mesh_.nEdges()),
633  nChangedEdges_(0),
634  nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
635  nEvals_(0),
636  nUnvisitedPoints_(mesh_.nPoints()),
637  nUnvisitedEdges_(mesh_.nEdges())
638 {
639  if (allPointInfo_.size() != mesh_.nPoints())
640  {
642  (
643  "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
644  "(const polyMesh&, const labelList&, const List<Type>,"
645  " List<Type>&, List<Type>&, const label maxIter)"
646  ) << "size of pointInfo work array is not equal to the number"
647  << " of points in the mesh" << endl
648  << " pointInfo :" << allPointInfo_.size() << endl
649  << " mesh.nPoints:" << mesh_.nPoints()
650  << exit(FatalError);
651  }
652  if (allEdgeInfo_.size() != mesh_.nEdges())
653  {
655  (
656  "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
657  "(const polyMesh&, const labelList&, const List<Type>,"
658  " List<Type>&, List<Type>&, const label maxIter)"
659  ) << "size of edgeInfo work array is not equal to the number"
660  << " of edges in the mesh" << endl
661  << " edgeInfo :" << allEdgeInfo_.size() << endl
662  << " mesh.nEdges:" << mesh_.nEdges()
663  << exit(FatalError);
664  }
665 
666 
667  // Set from initial changed points data
668  setPointInfo(changedPoints, changedPointsInfo);
669 
670  if (debug)
671  {
672  Info<< typeName << ": Seed points : "
673  << returnReduce(nChangedPoints_, sumOp<label>()) << endl;
674  }
675 
676  // Iterate until nothing changes
677  label iter = iterate(maxIter);
678 
679  if ((maxIter > 0) && (iter >= maxIter))
680  {
682  (
683  "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
684  "(const polyMesh&, const labelList&, const List<Type>,"
685  " List<Type>&, List<Type>&, const label maxIter)"
686  ) << "Maximum number of iterations reached. Increase maxIter." << endl
687  << " maxIter:" << maxIter << endl
688  << " nChangedPoints:" << nChangedPoints_ << endl
689  << " nChangedEdges:" << nChangedEdges_ << endl
690  << exit(FatalError);
691  }
692 }
693 
694 
695 template<class Type, class TrackingData>
697 (
698  const polyMesh& mesh,
699  UList<Type>& allPointInfo,
700  UList<Type>& allEdgeInfo,
701  TrackingData& td
702 )
703 :
704  mesh_(mesh),
705  allPointInfo_(allPointInfo),
706  allEdgeInfo_(allEdgeInfo),
707  td_(td),
708  changedPoint_(mesh_.nPoints(), false),
709  changedPoints_(mesh_.nPoints()),
710  nChangedPoints_(0),
711  changedEdge_(mesh_.nEdges(), false),
712  changedEdges_(mesh_.nEdges()),
713  nChangedEdges_(0),
714  nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
715  nEvals_(0),
716  nUnvisitedPoints_(mesh_.nPoints()),
717  nUnvisitedEdges_(mesh_.nEdges())
718 {}
719 
720 
721 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
722 
723 template<class Type, class TrackingData>
725 {}
726 
727 
728 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
729 
730 
731 template<class Type, class TrackingData>
733 {
734  return nUnvisitedPoints_;
735 }
736 
737 
738 template<class Type, class TrackingData>
740 {
741  return nUnvisitedEdges_;
742 }
743 
744 
745 // Copy point information into member data
746 template<class Type, class TrackingData>
748 (
749  const labelList& changedPoints,
750  const List<Type>& changedPointsInfo
751 )
752 {
753  forAll(changedPoints, changedPointI)
754  {
755  label pointI = changedPoints[changedPointI];
756 
757  bool wasValid = allPointInfo_[pointI].valid(td_);
758 
759  // Copy info for pointI
760  allPointInfo_[pointI] = changedPointsInfo[changedPointI];
761 
762  // Maintain count of unset points
763  if (!wasValid && allPointInfo_[pointI].valid(td_))
764  {
765  --nUnvisitedPoints_;
766  }
767 
768  // Mark pointI as changed, both on list and on point itself.
769 
770  if (!changedPoint_[pointI])
771  {
772  changedPoint_[pointI] = true;
773  changedPoints_[nChangedPoints_++] = pointI;
774  }
775  }
776 
777  // Sync
778  handleCollocatedPoints();
779 }
780 
781 
782 // Propagate information from edge to point. Return number of points changed.
783 template<class Type, class TrackingData>
785 {
786  for
787  (
788  label changedEdgeI = 0;
789  changedEdgeI < nChangedEdges_;
790  changedEdgeI++
791  )
792  {
793  label edgeI = changedEdges_[changedEdgeI];
794 
795  if (!changedEdge_[edgeI])
796  {
797  FatalErrorIn("PointEdgeWave<Type, TrackingData>::edgeToPoint()")
798  << "edge " << edgeI
799  << " not marked as having been changed" << nl
800  << "This might be caused by multiple occurences of the same"
801  << " seed point." << abort(FatalError);
802  }
803 
804 
805  const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
806 
807  // Evaluate all connected points (= edge endpoints)
808  const edge& e = mesh_.edges()[edgeI];
809 
810  forAll(e, eI)
811  {
812  Type& currentWallInfo = allPointInfo_[e[eI]];
813 
814  if (!currentWallInfo.equal(neighbourWallInfo, td_))
815  {
816  updatePoint
817  (
818  e[eI],
819  edgeI,
820  neighbourWallInfo,
821  currentWallInfo
822  );
823  }
824  }
825 
826  // Reset status of edge
827  changedEdge_[edgeI] = false;
828  }
829 
830  // Handled all changed edges by now
831  nChangedEdges_ = 0;
832 
833  if (nCyclicPatches_ > 0)
834  {
835  // Transfer changed points across cyclic halves
836  handleCyclicPatches();
837  }
838  if (Pstream::parRun())
839  {
840  // Transfer changed points from neighbouring processors.
841  handleProcPatches();
842  }
843 
844  //if (debug)
845  //{
846  // Pout<< "Changed points : " << nChangedPoints_ << endl;
847  //}
848 
849  // Sum nChangedPoints over all procs
850  label totNChanged = nChangedPoints_;
851 
852  reduce(totNChanged, sumOp<label>());
853 
854  return totNChanged;
855 }
856 
857 
858 // Propagate information from point to edge. Return number of edges changed.
859 template<class Type, class TrackingData>
861 {
862  const labelListList& pointEdges = mesh_.pointEdges();
863 
864  for
865  (
866  label changedPointI = 0;
867  changedPointI < nChangedPoints_;
868  changedPointI++
869  )
870  {
871  label pointI = changedPoints_[changedPointI];
872 
873  if (!changedPoint_[pointI])
874  {
875  FatalErrorIn("PointEdgeWave<Type, TrackingData>::pointToEdge()")
876  << "Point " << pointI
877  << " not marked as having been changed" << nl
878  << "This might be caused by multiple occurences of the same"
879  << " seed point." << abort(FatalError);
880  }
881 
882  const Type& neighbourWallInfo = allPointInfo_[pointI];
883 
884  // Evaluate all connected edges
885 
886  const labelList& edgeLabels = pointEdges[pointI];
887  forAll(edgeLabels, edgeLabelI)
888  {
889  label edgeI = edgeLabels[edgeLabelI];
890 
891  Type& currentWallInfo = allEdgeInfo_[edgeI];
892 
893  if (!currentWallInfo.equal(neighbourWallInfo, td_))
894  {
895  updateEdge
896  (
897  edgeI,
898  pointI,
899  neighbourWallInfo,
900  currentWallInfo
901  );
902  }
903  }
904 
905  // Reset status of point
906  changedPoint_[pointI] = false;
907  }
908 
909  // Handled all changed points by now
910  nChangedPoints_ = 0;
911 
912  //if (debug)
913  //{
914  // Pout<< "Changed edges : " << nChangedEdges_ << endl;
915  //}
916 
917  // Sum nChangedPoints over all procs
918  label totNChanged = nChangedEdges_;
919 
920  reduce(totNChanged, sumOp<label>());
921 
922  return totNChanged;
923 }
924 
925 
926 // Iterate
927 template<class Type, class TrackingData>
929 (
930  const label maxIter
931 )
932 {
933  if (nCyclicPatches_ > 0)
934  {
935  // Transfer changed points across cyclic halves
936  handleCyclicPatches();
937  }
938  if (Pstream::parRun())
939  {
940  // Transfer changed points from neighbouring processors.
941  handleProcPatches();
942  }
943 
944  nEvals_ = 0;
945 
946  label iter = 0;
947 
948  while (iter < maxIter)
949  {
950  while (iter < maxIter)
951  {
952  if (debug)
953  {
954  Info<< typeName << ": Iteration " << iter << endl;
955  }
956 
957  label nEdges = pointToEdge();
958 
959  if (debug)
960  {
961  Info<< typeName << ": Total changed edges : "
962  << nEdges << endl;
963  }
964 
965  if (nEdges == 0)
966  {
967  break;
968  }
969 
970  label nPoints = edgeToPoint();
971 
972  if (debug)
973  {
974  Info<< typeName << ": Total changed points : "
975  << nPoints << nl
976  << typeName << ": Total evaluations : "
977  << returnReduce(nEvals_, sumOp<label>()) << nl
978  << typeName << ": Remaining unvisited points: "
979  << returnReduce(nUnvisitedPoints_, sumOp<label>()) << nl
980  << typeName << ": Remaining unvisited edges : "
981  << returnReduce(nUnvisitedEdges_, sumOp<label>()) << nl
982  << endl;
983  }
984 
985  if (nPoints == 0)
986  {
987  break;
988  }
989 
990  iter++;
991  }
992 
993 
994  // Enforce collocated points are exactly equal. This might still mean
995  // non-collocated points are not equal though. WIP.
996  label nPoints = handleCollocatedPoints();
997  if (debug)
998  {
999  Info<< typeName << ": Collocated point sync : "
1000  << nPoints << nl << endl;
1001  }
1002 
1003  if (nPoints == 0)
1004  {
1005  break;
1006  }
1007  }
1008 
1009  return iter;
1010 }
1011 
1012 
1013 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:380
const labelListList & globalPointSlaves() const
virtual const tensorField & forwardT() const
Return face transformation tensor.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
combineEqOp(TrackingData &td)
Definition: PointEdgeWave.C:55
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
const word & name() const
Return name.
static void distribute(const Pstream::commsTypes commsType, const List< labelPair > &schedule, const label constructSize, const labelListList &subMap, const labelListList &constructMap, List< T > &, const int tag=UPstream::msgType())
Distribute data. Note:schedule only used for Pstream::scheduled.
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Class containing processor-to-processor mapping information.
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:242
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
label pointToEdge()
Propagate from point to edge. Returns total number of edges.
label getUnsetEdges() const
Get number of unvisited edges, i.e. edges that were not (yet)
messageStream Info
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:278
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
label nPoints() const
Return number of points supporting patch faces.
Namespace for OpenFOAM.
void setPointInfo(const labelList &changedPoints, const List< Type > &changedPointsInfo)
Copy initial data into allPointInfo_.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label getUnsetPoints() const
void reverseDistribute(const label constructSize, List< T > &, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Reverse distribute data using default commsType.
static const char nl
Definition: Ostream.H:260
const double e
Elementary charge.
Definition: doubleFloat.H:78
const mapDistribute & globalPointSlavesMap() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Neighbour processor patch.
virtual bool parallel() const
Are the cyclic planes parallel.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
#define forAll(list, i)
Definition: UList.H:421
label edgeToPoint()
Propagate from edge to point. Returns total number of points.
Cyclic plane patch.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label iterate(const label maxIter)
Iterate until no changes or maxIter reached. Returns actual.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A list of faces which address into the list of points.
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
error FatalError
scalar y
label nPatches
Definition: readKivaGrid.H:402
void operator()(Type &x, const Type &y) const
Definition: PointEdgeWave.C:60
const cyclicPolyPatch & neighbPatch() const
Wave propagation of information through grid. Every iteration information goes through one layer of e...
Definition: PointEdgeWave.H:85
int neighbProcNo() const
Return neigbour processor number.
const edgeList & coupledPoints() const
Return connected points (from patch local to neighbour patch local)
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
label nPoints
Reduction class. If x and y are not equal assign value.
Definition: PointEdgeWave.C:50
const labelList & neighbPoints() const
Return neighbour point labels. WIP.
void reserve(const label)
Reserve allocation space for at least this size.
Definition: DynamicListI.H:130
bool equal(const T &s1, const T &s2)
Definition: doubleFloat.H:62
~PointEdgeWave()
Destructor.
const Field< PointType > & points() const
Return reference to global points.
Combination-Reduction operation for a parallel run. The information from all nodes is collected on th...
label constructSize() const
Constructed data size.