FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "PointEdgeWave.H"
27 #include <OpenFOAM/polyMesh.H>
30 #include <OpenFOAM/OPstream.H>
31 #include <OpenFOAM/IPstream.H>
33 #include <OpenFOAM/debug.H>
34 #include <OpenFOAM/typeInfo.H>
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 template <class Type>
40 
41 
42 // Offset labelList. Used for transferring from one cyclic half to the other.
43 template <class Type>
44 void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
45 {
46  forAll(elems, i)
47  {
48  elems[i] += val;
49  }
50 }
51 
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
55 // Gets point-point correspondence. Is
56 // - list of halfA points (in cyclic patch points)
57 // - list of halfB points (can overlap with A!)
58 // - for every patchPoint its corresponding point
59 template <class Type>
61 {
62  label cycHalf = 0;
63 
64  forAll(mesh_.boundaryMesh(), patchI)
65  {
66  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
67 
68  if (isA<cyclicPolyPatch>(patch))
69  {
70  label halfSize = patch.size()/2;
71 
72  SubList<face> halfAFaces
73  (
74  mesh_.faces(),
75  halfSize,
76  patch.start()
77  );
78 
79  cycHalves_.set
80  (
81  cycHalf++,
82  new primitivePatch(halfAFaces, mesh_.points())
83  );
84 
85  SubList<face> halfBFaces
86  (
87  mesh_.faces(),
88  halfSize,
89  patch.start() + halfSize
90  );
91 
92  cycHalves_.set
93  (
94  cycHalf++,
95  new primitivePatch(halfBFaces, mesh_.points())
96  );
97  }
98  }
99 }
100 
101 
102 // Handle leaving domain. Implementation referred to Type
103 template <class Type>
105 (
106  const polyPatch& meshPatch,
107  const primitivePatch& patch,
108  const labelList& patchPointLabels,
109  List<Type>& pointInfo
110 ) const
111 {
112  const labelList& meshPoints = patch.meshPoints();
113 
114  forAll(patchPointLabels, i)
115  {
116  label patchPointI = patchPointLabels[i];
117 
118  const point& pt = patch.points()[meshPoints[patchPointI]];
119 
120  pointInfo[i].leaveDomain(meshPatch, patchPointI, pt);
121  }
122 }
123 
124 
125 // Handle entering domain. Implementation referred to Type
126 template <class Type>
128 (
129  const polyPatch& meshPatch,
130  const primitivePatch& patch,
131  const labelList& patchPointLabels,
132  List<Type>& pointInfo
133 ) const
134 {
135  const labelList& meshPoints = patch.meshPoints();
136 
137  forAll(patchPointLabels, i)
138  {
139  label patchPointI = patchPointLabels[i];
140 
141  const point& pt = patch.points()[meshPoints[patchPointI]];
142 
143  pointInfo[i].enterDomain(meshPatch, patchPointI, pt);
144  }
145 }
146 
147 
148 // Transform. Implementation referred to Type
149 template <class Type>
151 (
152  const tensorField& rotTensor,
153  List<Type>& pointInfo
154 ) const
155 {
156  if (rotTensor.size() == 1)
157  {
158  const tensor& T = rotTensor[0];
159 
160  forAll(pointInfo, i)
161  {
162  pointInfo[i].transform(T);
163  }
164  }
165  else
166  {
168  (
169  "PointEdgeWave<Type>::transform(const tensorField&, List<Type>&)"
170  ) << "Parallel cyclics not supported" << abort(FatalError);
171 
172  forAll(pointInfo, i)
173  {
174  pointInfo[i].transform(rotTensor[i]);
175  }
176  }
177 }
178 
179 
180 // Update info for pointI, at position pt, with information from
181 // neighbouring edge.
182 // Updates:
183 // - changedPoint_, changedPoints_, nChangedPoints_,
184 // - statistics: nEvals_, nUnvisitedPoints_
185 template <class Type>
187 (
188  const label pointI,
189  const label neighbourEdgeI,
190  const Type& neighbourInfo,
191  const scalar tol,
192  Type& pointInfo
193 )
194 {
195  nEvals_++;
196 
197  bool wasValid = pointInfo.valid();
198 
199  bool propagate =
200  pointInfo.updatePoint
201  (
202  mesh_,
203  pointI,
204  neighbourEdgeI,
205  neighbourInfo,
206  tol
207  );
208 
209  if (propagate)
210  {
211  if (!changedPoint_[pointI])
212  {
213  changedPoint_[pointI] = true;
214  changedPoints_[nChangedPoints_++] = pointI;
215  }
216  }
217 
218  if (!wasValid && pointInfo.valid())
219  {
220  --nUnvisitedPoints_;
221  }
222 
223  return propagate;
224 }
225 
226 
227 // Update info for pointI, at position pt, with information from
228 // same point.
229 // Updates:
230 // - changedPoint_, changedPoints_, nChangedPoints_,
231 // - statistics: nEvals_, nUnvisitedPoints_
232 template <class Type>
234 (
235  const label pointI,
236  const Type& neighbourInfo,
237  const scalar tol,
238  Type& pointInfo
239 )
240 {
241  nEvals_++;
242 
243  bool wasValid = pointInfo.valid();
244 
245  bool propagate =
246  pointInfo.updatePoint
247  (
248  mesh_,
249  pointI,
250  neighbourInfo,
251  tol
252  );
253 
254  if (propagate)
255  {
256  if (!changedPoint_[pointI])
257  {
258  changedPoint_[pointI] = true;
259  changedPoints_[nChangedPoints_++] = pointI;
260  }
261  }
262 
263  if (!wasValid && pointInfo.valid())
264  {
265  --nUnvisitedPoints_;
266  }
267 
268  return propagate;
269 }
270 
271 
272 // Update info for edgeI, at position pt, with information from
273 // neighbouring point.
274 // Updates:
275 // - changedEdge_, changedEdges_, nChangedEdges_,
276 // - statistics: nEvals_, nUnvisitedEdge_
277 template <class Type>
279 (
280  const label edgeI,
281  const label neighbourPointI,
282  const Type& neighbourInfo,
283  const scalar tol,
284  Type& edgeInfo
285 )
286 {
287  nEvals_++;
288 
289  bool wasValid = edgeInfo.valid();
290 
291  bool propagate =
292  edgeInfo.updateEdge
293  (
294  mesh_,
295  edgeI,
296  neighbourPointI,
297  neighbourInfo,
298  tol
299  );
300 
301  if (propagate)
302  {
303  if (!changedEdge_[edgeI])
304  {
305  changedEdge_[edgeI] = true;
306  changedEdges_[nChangedEdges_++] = edgeI;
307  }
308  }
309 
310  if (!wasValid && edgeInfo.valid())
311  {
312  --nUnvisitedEdges_;
313  }
314 
315  return propagate;
316 }
317 
318 
319 // Check if patches of given type name are present
320 template <class Type>
321 template <class PatchType>
323 {
324  label nPatches = 0;
325 
326  forAll(mesh_.boundaryMesh(), patchI)
327  {
328  if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
329  {
330  nPatches++;
331  }
332  }
333  return nPatches;
334 }
335 
336 
337 // Collect changed patch points
338 template <class Type>
340 (
341  const primitivePatch& patch,
342 
343  DynamicList<Type>& patchInfo,
344  DynamicList<label>& patchPoints,
345  DynamicList<label>& owner,
346  DynamicList<label>& ownerIndex
347 ) const
348 {
349  const labelList& meshPoints = patch.meshPoints();
350  const faceList& localFaces = patch.localFaces();
351  const labelListList& pointFaces = patch.pointFaces();
352 
353  forAll(meshPoints, patchPointI)
354  {
355  label meshPointI = meshPoints[patchPointI];
356 
357  if (changedPoint_[meshPointI])
358  {
359  patchInfo.append(allPointInfo_[meshPointI]);
360  patchPoints.append(patchPointI);
361 
362  label patchFaceI = pointFaces[patchPointI][0];
363 
364  const face& f = localFaces[patchFaceI];
365 
366  label index = findIndex(f, patchPointI);
367 
368  owner.append(patchFaceI);
369  ownerIndex.append(index);
370  }
371  }
372 
373  patchInfo.shrink();
374  patchPoints.shrink();
375  owner.shrink();
376  ownerIndex.shrink();
377 }
378 
379 
380 // Update overall for changed patch points
381 template <class Type>
383 (
384  const polyPatch& meshPatch,
385  const primitivePatch& patch,
386  const labelList& owner,
387  const labelList& ownerIndex,
388  List<Type>& patchInfo
389 )
390 {
391  const faceList& localFaces = patch.localFaces();
392  const labelList& meshPoints = patch.meshPoints();
393 
394  // Get patch and mesh points.
395  labelList changedPatchPoints(patchInfo.size());
396  labelList changedMeshPoints(patchInfo.size());
397 
398  forAll(owner, i)
399  {
400  label faceI = owner[i];
401 
402  const face& f = localFaces[faceI];
403 
404  label index = (f.size() - ownerIndex[i]) % f.size();
405 
406  changedPatchPoints[i] = f[index];
407  changedMeshPoints[i] = meshPoints[f[index]];
408  }
409 
410  // Adapt for entering domain
411  enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
412 
413  // Merge received info
414  forAll(patchInfo, i)
415  {
416  updatePoint
417  (
418  changedMeshPoints[i],
419  patchInfo[i],
420  propagationTol_,
421  allPointInfo_[changedMeshPoints[i]]
422  );
423  }
424 }
425 
426 
427 
428 // Transfer all the information to/from neighbouring processors
429 template <class Type>
431 {
432  // 1. Send all point info on processor patches. Send as
433  // face label + offset in face.
434 
435  forAll(mesh_.boundaryMesh(), patchI)
436  {
437  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
438 
439  if (isA<processorPolyPatch>(patch))
440  {
441  // Get all changed points in relative addressing
442 
443  DynamicList<Type> patchInfo(patch.nPoints());
444  DynamicList<label> patchPoints(patch.nPoints());
445  DynamicList<label> owner(patch.nPoints());
446  DynamicList<label> ownerIndex(patch.nPoints());
447 
448  getChangedPatchPoints
449  (
450  patch,
451  patchInfo,
452  patchPoints,
453  owner,
454  ownerIndex
455  );
456 
457  // Adapt for leaving domain
458  leaveDomain(patch, patch, patchPoints, patchInfo);
459 
460  const processorPolyPatch& procPatch =
461  refCast<const processorPolyPatch>(patch);
462 
463  if (debug)
464  {
465  Pout<< "Processor patch " << patchI << ' ' << patch.name()
466  << " communicating with " << procPatch.neighbProcNo()
467  << " Sending:" << patchInfo.size() << endl;
468  }
469 
470  {
471  OPstream toNeighbour
472  (
473  Pstream::blocking,
474  procPatch.neighbProcNo()
475  );
476 
477  toNeighbour << owner << ownerIndex << patchInfo;
478  }
479  }
480  }
481 
482 
483  //
484  // 2. Receive all point info on processor patches.
485  //
486 
487  forAll(mesh_.boundaryMesh(), patchI)
488  {
489  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
490 
491  if (isA<processorPolyPatch>(patch))
492  {
493  const processorPolyPatch& procPatch =
494  refCast<const processorPolyPatch>(patch);
495 
496  List<Type> patchInfo;
497  labelList owner;
498  labelList ownerIndex;
499  {
500  IPstream fromNeighbour
501  (
502  Pstream::blocking,
503  procPatch.neighbProcNo()
504  );
505 
506  fromNeighbour >> owner >> ownerIndex >> patchInfo;
507  }
508 
509  if (debug)
510  {
511  Pout<< "Processor patch " << patchI << ' ' << patch.name()
512  << " communicating with " << procPatch.neighbProcNo()
513  << " Received:" << patchInfo.size() << endl;
514  }
515 
516  // Apply transform to received data for non-parallel planes
517  if (!procPatch.parallel())
518  {
519  transform(procPatch.forwardT(), patchInfo);
520  }
521 
522  updateFromPatchInfo
523  (
524  patch,
525  patch,
526  owner,
527  ownerIndex,
528  patchInfo
529  );
530  }
531  }
532 
533 
534 
535  //
536  // 3. Handle all shared points
537  // (Note:irrespective if changed or not for now)
538  //
539 
540  const globalMeshData& pd = mesh_.globalData();
541 
542  List<Type> sharedData(pd.nGlobalPoints());
543 
544  forAll(pd.sharedPointLabels(), i)
545  {
546  label meshPointI = pd.sharedPointLabels()[i];
547 
548  // Fill my entries in the shared points
549  sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
550  }
551 
552  // Combine on master. Reduce operator has to handle a list and call
553  // Type.updatePoint for all elements
554  combineReduce(sharedData, listUpdateOp<Type>());
555 
556  forAll(pd.sharedPointLabels(), i)
557  {
558  label meshPointI = pd.sharedPointLabels()[i];
559 
560  // Retrieve my entries from the shared points
561  updatePoint
562  (
563  meshPointI,
564  sharedData[pd.sharedPointAddr()[i]],
565  propagationTol_,
566  allPointInfo_[meshPointI]
567  );
568  }
569 }
570 
571 
572 template <class Type>
574 {
575  // 1. Send all point info on cyclic patches. Send as
576  // face label + offset in face.
577 
578  label cycHalf = 0;
579 
580  forAll(mesh_.boundaryMesh(), patchI)
581  {
582  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
583 
584  if (isA<cyclicPolyPatch>(patch))
585  {
586  const primitivePatch& halfA = cycHalves_[cycHalf++];
587  const primitivePatch& halfB = cycHalves_[cycHalf++];
588 
589  // HalfA : get all changed points in relative addressing
590 
591  DynamicList<Type> halfAInfo(halfA.nPoints());
592  DynamicList<label> halfAPoints(halfA.nPoints());
593  DynamicList<label> halfAOwner(halfA.nPoints());
594  DynamicList<label> halfAIndex(halfA.nPoints());
595 
596  getChangedPatchPoints
597  (
598  halfA,
599  halfAInfo,
600  halfAPoints,
601  halfAOwner,
602  halfAIndex
603  );
604 
605  // HalfB : get all changed points in relative addressing
606 
607  DynamicList<Type> halfBInfo(halfB.nPoints());
608  DynamicList<label> halfBPoints(halfB.nPoints());
609  DynamicList<label> halfBOwner(halfB.nPoints());
610  DynamicList<label> halfBIndex(halfB.nPoints());
611 
612  getChangedPatchPoints
613  (
614  halfB,
615  halfBInfo,
616  halfBPoints,
617  halfBOwner,
618  halfBIndex
619  );
620 
621 
622  // HalfA : adapt for leaving domain
623  leaveDomain(patch, halfA, halfAPoints, halfAInfo);
624 
625  // HalfB : adapt for leaving domain
626  leaveDomain(patch, halfB, halfBPoints, halfBInfo);
627 
628 
629  // Apply rotation for non-parallel planes
630  const cyclicPolyPatch& cycPatch =
631  refCast<const cyclicPolyPatch>(patch);
632 
633  if (!cycPatch.parallel())
634  {
635  // received data from half1
636  transform(cycPatch.reverseT(), halfAInfo);
637 
638  // received data from half2
639  transform(cycPatch.forwardT(), halfBInfo);
640  }
641 
642  if (debug)
643  {
644  Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
645  << " Changed on first half : " << halfAInfo.size()
646  << " Changed on second half : " << halfBInfo.size()
647  << endl;
648  }
649 
650  // Half1: update with data from halfB
651  updateFromPatchInfo
652  (
653  patch,
654  halfA,
655  halfBOwner,
656  halfBIndex,
657  halfBInfo
658  );
659 
660  // Half2: update with data from halfA
661  updateFromPatchInfo
662  (
663  patch,
664  halfB,
665  halfAOwner,
666  halfAIndex,
667  halfAInfo
668  );
669 
670  if (debug)
671  {
672  //checkCyclic(patch);
673  }
674  }
675  }
676 }
677 
678 
679 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
680 
681 // Iterate, propagating changedPointsInfo across mesh, until no change (or
682 // maxIter reached). Initial point values specified.
683 template <class Type>
685 (
686  const polyMesh& mesh,
687  const labelList& changedPoints,
688  const List<Type>& changedPointsInfo,
689 
690  List<Type>& allPointInfo,
691  List<Type>& allEdgeInfo,
692 
693  const label maxIter
694 )
695 :
696  mesh_(mesh),
697  allPointInfo_(allPointInfo),
698  allEdgeInfo_(allEdgeInfo),
699  changedPoint_(mesh_.nPoints(), false),
700  changedPoints_(mesh_.nPoints()),
701  nChangedPoints_(0),
702  changedEdge_(mesh_.nEdges(), false),
703  changedEdges_(mesh_.nEdges()),
704  nChangedEdges_(0),
705  nCyclicPatches_(countPatchType<cyclicPolyPatch>()),
706  cycHalves_(2*nCyclicPatches_),
707  nEvals_(0),
708  nUnvisitedPoints_(mesh_.nPoints()),
709  nUnvisitedEdges_(mesh_.nEdges())
710 {
711  if (allPointInfo_.size() != mesh_.nPoints())
712  {
714  (
715  "PointEdgeWave<Type>::PointEdgeWave"
716  "(const polyMesh&, const labelList&, const List<Type>,"
717  " List<Type>&, List<Type>&, const label maxIter)"
718  ) << "size of pointInfo work array is not equal to the number"
719  << " of points in the mesh" << endl
720  << " pointInfo :" << allPointInfo_.size() << endl
721  << " mesh.nPoints:" << mesh_.nPoints()
722  << exit(FatalError);
723  }
724  if (allEdgeInfo_.size() != mesh_.nEdges())
725  {
727  (
728  "PointEdgeWave<Type>::PointEdgeWave"
729  "(const polyMesh&, const labelList&, const List<Type>,"
730  " List<Type>&, List<Type>&, const label maxIter)"
731  ) << "size of edgeInfo work array is not equal to the number"
732  << " of edges in the mesh" << endl
733  << " edgeInfo :" << allEdgeInfo_.size() << endl
734  << " mesh.nEdges:" << mesh_.nEdges()
735  << exit(FatalError);
736  }
737 
738 
739  // Calculate cyclic halves addressing.
740  if (nCyclicPatches_ > 0)
741  {
742  calcCyclicAddressing();
743  }
744 
745 
746  // Set from initial changed points data
747  setPointInfo(changedPoints, changedPointsInfo);
748 
749  if (debug)
750  {
751  Pout<< "Seed points : " << nChangedPoints_ << endl;
752  }
753 
754  // Iterate until nothing changes
755  label iter = iterate(maxIter);
756 
757  if ((maxIter > 0) && (iter >= maxIter))
758  {
760  (
761  "PointEdgeWave<Type>::PointEdgeWave"
762  "(const polyMesh&, const labelList&, const List<Type>,"
763  " List<Type>&, List<Type>&, const label maxIter)"
764  ) << "Maximum number of iterations reached. Increase maxIter." << endl
765  << " maxIter:" << maxIter << endl
766  << " nChangedPoints:" << nChangedPoints_ << endl
767  << " nChangedEdges:" << nChangedEdges_ << endl
768  << exit(FatalError);
769  }
770 }
771 
772 
773 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
774 
775 template <class Type>
777 {}
778 
779 
780 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
781 
782 
783 template <class Type>
785 {
786  return nUnvisitedPoints_;
787 }
788 
789 
790 template <class Type>
792 {
793  return nUnvisitedEdges_;
794 }
795 
796 
797 // Copy point information into member data
798 template <class Type>
800 (
801  const labelList& changedPoints,
802  const List<Type>& changedPointsInfo
803 )
804 {
805  forAll(changedPoints, changedPointI)
806  {
807  label pointI = changedPoints[changedPointI];
808 
809  bool wasValid = allPointInfo_[pointI].valid();
810 
811  // Copy info for pointI
812  allPointInfo_[pointI] = changedPointsInfo[changedPointI];
813 
814  // Maintain count of unset points
815  if (!wasValid && allPointInfo_[pointI].valid())
816  {
817  --nUnvisitedPoints_;
818  }
819 
820  // Mark pointI as changed, both on list and on point itself.
821 
822  if (!changedPoint_[pointI])
823  {
824  changedPoint_[pointI] = true;
825  changedPoints_[nChangedPoints_++] = pointI;
826  }
827  }
828 }
829 
830 
831 // Propagate information from edge to point. Return number of points changed.
832 template <class Type>
834 {
835  for
836  (
837  label changedEdgeI = 0;
838  changedEdgeI < nChangedEdges_;
839  changedEdgeI++
840  )
841  {
842  label edgeI = changedEdges_[changedEdgeI];
843 
844  if (!changedEdge_[edgeI])
845  {
846  FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
847  << "edge " << edgeI
848  << " not marked as having been changed" << nl
849  << "This might be caused by multiple occurences of the same"
850  << " seed point." << abort(FatalError);
851  }
852 
853 
854  const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
855 
856  // Evaluate all connected points (= edge endpoints)
857  const edge& e = mesh_.edges()[edgeI];
858 
859  forAll(e, eI)
860  {
861  Type& currentWallInfo = allPointInfo_[e[eI]];
862 
863  if (currentWallInfo != neighbourWallInfo)
864  {
865  updatePoint
866  (
867  e[eI],
868  edgeI,
869  neighbourWallInfo,
870  propagationTol_,
871  currentWallInfo
872  );
873  }
874  }
875 
876  // Reset status of edge
877  changedEdge_[edgeI] = false;
878  }
879 
880  // Handled all changed edges by now
881  nChangedEdges_ = 0;
882 
883  if (nCyclicPatches_ > 0)
884  {
885  // Transfer changed points across cyclic halves
886  handleCyclicPatches();
887  }
888  if (Pstream::parRun())
889  {
890  // Transfer changed points from neighbouring processors.
891  handleProcPatches();
892  }
893 
894  if (debug)
895  {
896  Pout<< "Changed points : " << nChangedPoints_ << endl;
897  }
898 
899  // Sum nChangedPoints over all procs
900  label totNChanged = nChangedPoints_;
901 
902  reduce(totNChanged, sumOp<label>());
903 
904  return totNChanged;
905 }
906 
907 
908 // Propagate information from point to edge. Return number of edges changed.
909 template <class Type>
911 {
912  const labelListList& pointEdges = mesh_.pointEdges();
913 
914  for
915  (
916  label changedPointI = 0;
917  changedPointI < nChangedPoints_;
918  changedPointI++
919  )
920  {
921  label pointI = changedPoints_[changedPointI];
922 
923  if (!changedPoint_[pointI])
924  {
925  FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
926  << "Point " << pointI
927  << " not marked as having been changed" << nl
928  << "This might be caused by multiple occurences of the same"
929  << " seed point." << abort(FatalError);
930  }
931 
932  const Type& neighbourWallInfo = allPointInfo_[pointI];
933 
934  // Evaluate all connected edges
935 
936  const labelList& edgeLabels = pointEdges[pointI];
937  forAll(edgeLabels, edgeLabelI)
938  {
939  label edgeI = edgeLabels[edgeLabelI];
940 
941  Type& currentWallInfo = allEdgeInfo_[edgeI];
942 
943  if (currentWallInfo != neighbourWallInfo)
944  {
945  updateEdge
946  (
947  edgeI,
948  pointI,
949  neighbourWallInfo,
950  propagationTol_,
951  currentWallInfo
952  );
953  }
954  }
955 
956  // Reset status of point
957  changedPoint_[pointI] = false;
958  }
959 
960  // Handled all changed points by now
961  nChangedPoints_ = 0;
962 
963  if (debug)
964  {
965  Pout<< "Changed edges : " << nChangedEdges_ << endl;
966  }
967 
968  // Sum nChangedPoints over all procs
969  label totNChanged = nChangedEdges_;
970 
971  reduce(totNChanged, sumOp<label>());
972 
973  return totNChanged;
974 }
975 
976 
977 // Iterate
978 template <class Type>
979 Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
980 {
981  if (nCyclicPatches_ > 0)
982  {
983  // Transfer changed points across cyclic halves
984  handleCyclicPatches();
985  }
986  if (Pstream::parRun())
987  {
988  // Transfer changed points from neighbouring processors.
989  handleProcPatches();
990  }
991 
992  nEvals_ = 0;
993 
994  label iter = 0;
995 
996  while(iter < maxIter)
997  {
998  if (debug)
999  {
1000  Pout<< "Iteration " << iter << endl;
1001  }
1002 
1003  label nEdges = pointToEdge();
1004 
1005  if (debug)
1006  {
1007  Pout<< "Total changed edges : " << nEdges << endl;
1008  }
1009 
1010  if (nEdges == 0)
1011  {
1012  break;
1013  }
1014 
1015  label nPoints = edgeToPoint();
1016 
1017  if (debug)
1018  {
1019  Pout<< "Total changed points : " << nPoints << endl;
1020 
1021  Pout<< "Total evaluations : " << nEvals_ << endl;
1022 
1023  Pout<< "Remaining unvisited points: " << nUnvisitedPoints_ << endl;
1024 
1025  Pout<< "Remaining unvisited edges : " << nUnvisitedEdges_ << endl;
1026 
1027  Pout<< endl;
1028  }
1029 
1030  if (nPoints == 0)
1031  {
1032  break;
1033  }
1034 
1035  iter++;
1036  }
1037 
1038  return iter;
1039 }
1040 
1041 // ************************ vim: set sw=4 sts=4 et: ************************ //