FreeFOAM The Cross-Platform CFD Toolkit
FaceCellWave.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 "FaceCellWave.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>
39 const Foam::scalar Foam::FaceCellWave<Type>::geomTol_ = 1e-6;
40 
41 template <class Type>
43 
44 // Write to ostream
45 template <class Type>
47 (
48  const label nFaces,
49  const labelList& faceLabels,
50  const List<Type>& faceInfo,
51  Ostream& os
52 )
53 {
54  // Write list contents depending on data format
55  if (os.format() == IOstream::ASCII)
56  {
57  os << nFaces;
58 
59  for(label i = 0; i < nFaces; i++)
60  {
61  os << ' ' << faceLabels[i];
62  }
63  for(label i = 0; i < nFaces; i++)
64  {
65  os << ' ' << faceInfo[i];
66  }
67  }
68  else
69  {
70  os << nFaces;
71 
72  for(label i = 0; i < nFaces; i++)
73  {
74  os << faceLabels[i];
75  }
76  for(label i = 0; i < nFaces; i++)
77  {
78  os << faceInfo[i];
79  }
80  }
81  return os;
82 }
83 
84 
85 // Read from istream
86 template <class Type>
88 (
89  label& nFaces,
90  labelList& faceLabels,
91  List<Type>& faceInfo,
92  Istream& is
93 )
94 {
95  is >> nFaces;
96 
97  for(label i = 0; i < nFaces; i++)
98  {
99  is >> faceLabels[i];
100  }
101  for(label i = 0; i < nFaces; i++)
102  {
103  is >> faceInfo[i];
104  }
105  return is;
106 }
107 
108 
109 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
110 
111 // Update info for cellI, at position pt, with information from
112 // neighbouring face/cell.
113 // Updates:
114 // - changedCell_, changedCells_, nChangedCells_,
115 // - statistics: nEvals_, nUnvisitedCells_
116 template <class Type>
118 (
119  const label cellI,
120  const label neighbourFaceI,
121  const Type& neighbourInfo,
122  const scalar tol,
123  Type& cellInfo
124 )
125 {
126  nEvals_++;
127 
128  bool wasValid = cellInfo.valid();
129 
130  bool propagate =
131  cellInfo.updateCell
132  (
133  mesh_,
134  cellI,
135  neighbourFaceI,
136  neighbourInfo,
137  tol
138  );
139 
140  if (propagate)
141  {
142  if (!changedCell_[cellI])
143  {
144  changedCell_[cellI] = true;
145  changedCells_[nChangedCells_++] = cellI;
146  }
147  }
148 
149  if (!wasValid && cellInfo.valid())
150  {
151  --nUnvisitedCells_;
152  }
153 
154  return propagate;
155 }
156 
157 
158 // Update info for faceI, at position pt, with information from
159 // neighbouring face/cell.
160 // Updates:
161 // - changedFace_, changedFaces_, nChangedFaces_,
162 // - statistics: nEvals_, nUnvisitedFaces_
163 template <class Type>
165 (
166  const label faceI,
167  const label neighbourCellI,
168  const Type& neighbourInfo,
169  const scalar tol,
170  Type& faceInfo
171 )
172 {
173  nEvals_++;
174 
175  bool wasValid = faceInfo.valid();
176 
177  bool propagate =
178  faceInfo.updateFace
179  (
180  mesh_,
181  faceI,
182  neighbourCellI,
183  neighbourInfo,
184  tol
185  );
186 
187  if (propagate)
188  {
189  if (!changedFace_[faceI])
190  {
191  changedFace_[faceI] = true;
192  changedFaces_[nChangedFaces_++] = faceI;
193  }
194  }
195 
196  if (!wasValid && faceInfo.valid())
197  {
198  --nUnvisitedFaces_;
199  }
200 
201  return propagate;
202 }
203 
204 
205 // Update info for faceI, at position pt, with information from
206 // same face.
207 // Updates:
208 // - changedFace_, changedFaces_, nChangedFaces_,
209 // - statistics: nEvals_, nUnvisitedFaces_
210 template <class Type>
212 (
213  const label faceI,
214  const Type& neighbourInfo,
215  const scalar tol,
216  Type& faceInfo
217 )
218 {
219  nEvals_++;
220 
221  bool wasValid = faceInfo.valid();
222 
223  bool propagate =
224  faceInfo.updateFace
225  (
226  mesh_,
227  faceI,
228  neighbourInfo,
229  tol
230  );
231 
232  if (propagate)
233  {
234  if (!changedFace_[faceI])
235  {
236  changedFace_[faceI] = true;
237  changedFaces_[nChangedFaces_++] = faceI;
238  }
239  }
240 
241  if (!wasValid && faceInfo.valid())
242  {
243  --nUnvisitedFaces_;
244  }
245 
246  return propagate;
247 }
248 
249 
250 // For debugging: check status on both sides of cyclic
251 template <class Type>
252 void Foam::FaceCellWave<Type>::checkCyclic(const polyPatch& patch) const
253 {
254  label cycOffset = patch.size()/2;
255 
256  for(label patchFaceI = 0; patchFaceI < cycOffset; patchFaceI++)
257  {
258  label i1 = patch.start() + patchFaceI;
259  label i2 = i1 + cycOffset;
260 
261  if (!allFaceInfo_[i1].sameGeometry(mesh_, allFaceInfo_[i2], geomTol_))
262  {
263  FatalErrorIn("FaceCellWave<Type>::checkCyclic(const polyPatch&)")
264  << "problem: i:" << i1 << " otheri:" << i2
265  << " faceInfo:" << allFaceInfo_[i1]
266  << " otherfaceInfo:" << allFaceInfo_[i2]
267  << abort(FatalError);
268  }
269 
270  if (changedFace_[i1] != changedFace_[i2])
271  {
272  FatalErrorIn("FaceCellWave<Type>::checkCyclic(const polyPatch&)")
273  << " problem: i:" << i1 << " otheri:" << i2
274  << " faceInfo:" << allFaceInfo_[i1]
275  << " otherfaceInfo:" << allFaceInfo_[i2]
276  << " changedFace:" << changedFace_[i1]
277  << " otherchangedFace:" << changedFace_[i2]
278  << abort(FatalError);
279  }
280  }
281 }
282 
283 
284 // Check if has cyclic patches
285 template <class Type>
287 {
288  forAll(mesh_.boundaryMesh(), patchI)
289  {
290  if (isA<cyclicPolyPatch>(mesh_.boundaryMesh()[patchI]))
291  {
292  return true;
293  }
294  }
295  return false;
296 }
297 
298 
299 // Copy face information into member data
300 template <class Type>
302 (
303  const labelList& changedFaces,
304  const List<Type>& changedFacesInfo
305 )
306 {
307  forAll(changedFaces, changedFaceI)
308  {
309  label faceI = changedFaces[changedFaceI];
310 
311  bool wasValid = allFaceInfo_[faceI].valid();
312 
313  // Copy info for faceI
314  allFaceInfo_[faceI] = changedFacesInfo[changedFaceI];
315 
316  // Maintain count of unset faces
317  if (!wasValid && allFaceInfo_[faceI].valid())
318  {
319  --nUnvisitedFaces_;
320  }
321 
322  // Mark faceI as changed, both on list and on face itself.
323 
324  changedFace_[faceI] = true;
325  changedFaces_[nChangedFaces_++] = faceI;
326  }
327 }
328 
329 
330 // Merge face information into member data
331 template <class Type>
333 (
334  const polyPatch& patch,
335  const label nFaces,
336  const labelList& changedFaces,
337  const List<Type>& changedFacesInfo,
338  const bool
339 )
340 {
341  for(label changedFaceI = 0; changedFaceI < nFaces; changedFaceI++)
342  {
343  const Type& neighbourWallInfo = changedFacesInfo[changedFaceI];
344  label patchFaceI = changedFaces[changedFaceI];
345 
346  label meshFaceI = patch.start() + patchFaceI;
347 
348  Type& currentWallInfo = allFaceInfo_[meshFaceI];
349 
350  if (currentWallInfo != neighbourWallInfo)
351  {
352  updateFace
353  (
354  meshFaceI,
355  neighbourWallInfo,
356  propagationTol_,
357  currentWallInfo
358  );
359  }
360  }
361 }
362 
363 
364 // Construct compact patchFace change arrays for a (slice of a) single patch.
365 // changedPatchFaces in local patch numbering.
366 // Return length of arrays.
367 template <class Type>
369 (
370  const polyPatch& patch,
371  const label startFaceI,
372  const label nFaces,
373  labelList& changedPatchFaces,
374  List<Type>& changedPatchFacesInfo
375 ) const
376 {
377  label nChangedPatchFaces = 0;
378 
379  for(label i = 0; i < nFaces; i++)
380  {
381  label patchFaceI = i + startFaceI;
382 
383  label meshFaceI = patch.start() + patchFaceI;
384 
385  if (changedFace_[meshFaceI])
386  {
387  changedPatchFaces[nChangedPatchFaces] = patchFaceI;
388  changedPatchFacesInfo[nChangedPatchFaces] = allFaceInfo_[meshFaceI];
389  nChangedPatchFaces++;
390  }
391  }
392  return nChangedPatchFaces;
393 }
394 
395 
396 // Handle leaving domain. Implementation referred to Type
397 template <class Type>
399 (
400  const polyPatch& patch,
401  const label nFaces,
402  const labelList& faceLabels,
403  List<Type>& faceInfo
404 ) const
405 {
406  const vectorField& fc = mesh_.faceCentres();
407 
408  for(label i = 0; i < nFaces; i++)
409  {
410  label patchFaceI = faceLabels[i];
411 
412  label meshFaceI = patch.start() + patchFaceI;
413  faceInfo[i].leaveDomain(mesh_, patch, patchFaceI, fc[meshFaceI]);
414  }
415 }
416 
417 
418 // Handle entering domain. Implementation referred to Type
419 template <class Type>
421 (
422  const polyPatch& patch,
423  const label nFaces,
424  const labelList& faceLabels,
425  List<Type>& faceInfo
426 ) const
427 {
428  const vectorField& fc = mesh_.faceCentres();
429 
430  for(label i = 0; i < nFaces; i++)
431  {
432  label patchFaceI = faceLabels[i];
433 
434  label meshFaceI = patch.start() + patchFaceI;
435  faceInfo[i].enterDomain(mesh_, patch, patchFaceI, fc[meshFaceI]);
436  }
437 }
438 
439 
440 // Transform. Implementation referred to Type
441 template <class Type>
443 (
444  const tensorField& rotTensor,
445  const label nFaces,
446  List<Type>& faceInfo
447 )
448 {
449  if (rotTensor.size() == 1)
450  {
451  const tensor& T = rotTensor[0];
452 
453  for(label faceI = 0; faceI < nFaces; faceI++)
454  {
455  faceInfo[faceI].transform(mesh_, T);
456  }
457  }
458  else
459  {
460  for(label faceI = 0; faceI < nFaces; faceI++)
461  {
462  faceInfo[faceI].transform(mesh_, rotTensor[faceI]);
463  }
464  }
465 }
466 
467 
468 // Send face info to neighbour.
469 template <class Type>
471 (
472  const label neighbour,
473  const label nFaces,
474  const labelList& faceLabels,
475  const List<Type>& faceInfo
476 ) const
477 {
478  OPstream toNeighbour(Pstream::blocking, neighbour);
479 
480  writeFaces(nFaces, faceLabels, faceInfo, toNeighbour);
481 }
482 
483 
484 // Receive face info from neighbour
485 template <class Type>
487 (
488  const label neighbour,
489  labelList& faceLabels,
490  List<Type>& faceInfo
491 ) const
492 {
493  IPstream fromNeighbour(Pstream::blocking, neighbour);
494 
495  label nFaces = 0;
496  readFaces(nFaces, faceLabels, faceInfo, fromNeighbour);
497 
498  return nFaces;
499 }
500 
501 
502 // Offset mesh face. Used for transferring from one cyclic half to the other.
503 template <class Type>
505 (
506  const polyPatch&,
507  const label cycOffset,
508  const label nFaces,
509  labelList& faces
510 )
511 {
512  for(label faceI = 0; faceI < nFaces; faceI++)
513  {
514  faces[faceI] += cycOffset;
515  }
516 }
517 
518 
519 // Tranfer all the information to/from neighbouring processors
520 template <class Type>
522 {
523  // Send all
524 
525  forAll(mesh_.boundaryMesh(), patchI)
526  {
527  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
528 
529  if (isA<processorPolyPatch>(patch))
530  {
531  // Allocate buffers
532  label nSendFaces;
533  labelList sendFaces(patch.size());
534  List<Type> sendFacesInfo(patch.size());
535 
536  // Determine which faces changed on current patch
537  nSendFaces = getChangedPatchFaces
538  (
539  patch,
540  0,
541  patch.size(),
542  sendFaces,
543  sendFacesInfo
544  );
545 
546  // Adapt wallInfo for leaving domain
547  leaveDomain
548  (
549  patch,
550  nSendFaces,
551  sendFaces,
552  sendFacesInfo
553  );
554 
555  const processorPolyPatch& procPatch =
556  refCast<const processorPolyPatch>(patch);
557 
558  if (debug)
559  {
560  Pout<< " Processor patch " << patchI << ' ' << patch.name()
561  << " communicating with " << procPatch.neighbProcNo()
562  << " Sending:" << nSendFaces
563  << endl;
564  }
565 
566  sendPatchInfo
567  (
568  procPatch.neighbProcNo(),
569  nSendFaces,
570  sendFaces,
571  sendFacesInfo
572  );
573  }
574  }
575 
576  // Receive all
577 
578  forAll(mesh_.boundaryMesh(), patchI)
579  {
580  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
581 
582  if (isA<processorPolyPatch>(patch))
583  {
584  const processorPolyPatch& procPatch =
585  refCast<const processorPolyPatch>(patch);
586 
587  // Allocate buffers
588  label nReceiveFaces;
589  labelList receiveFaces(patch.size());
590  List<Type> receiveFacesInfo(patch.size());
591 
592  nReceiveFaces = receivePatchInfo
593  (
594  procPatch.neighbProcNo(),
595  receiveFaces,
596  receiveFacesInfo
597  );
598 
599  if (debug)
600  {
601  Pout<< " Processor patch " << patchI << ' ' << patch.name()
602  << " communicating with " << procPatch.neighbProcNo()
603  << " Receiving:" << nReceiveFaces
604  << endl;
605  }
606 
607  // Apply transform to received data for non-parallel planes
608  if (!procPatch.parallel())
609  {
610  transform
611  (
612  procPatch.forwardT(),
613  nReceiveFaces,
614  receiveFacesInfo
615  );
616  }
617 
618  // Adapt wallInfo for entering domain
619  enterDomain
620  (
621  patch,
622  nReceiveFaces,
623  receiveFaces,
624  receiveFacesInfo
625  );
626 
627  // Merge received info
628  mergeFaceInfo
629  (
630  patch,
631  nReceiveFaces,
632  receiveFaces,
633  receiveFacesInfo,
634  procPatch.parallel()
635  );
636  }
637  }
638 }
639 
640 
641 // Transfer information across cyclic halves.
642 // Note: Cyclic is two patches in one: one side from 0..size/2 and other
643 // side from size/2 .. size.
644 template <class Type>
646 {
647  forAll(mesh_.boundaryMesh(), patchI)
648  {
649  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
650 
651  if (isA<cyclicPolyPatch>(patch))
652  {
653  label halfSize = patch.size()/2;
654 
655  // Allocate buffers
656  label nSendFaces;
657  labelList sendFaces(halfSize);
658  List<Type> sendFacesInfo(halfSize);
659 
660  label nReceiveFaces;
661  labelList receiveFaces(halfSize);
662  List<Type> receiveFacesInfo(halfSize);
663 
664  // Half1: Determine which faces changed. Use sendFaces for storage
665  nSendFaces = getChangedPatchFaces
666  (
667  patch,
668  0,
669  halfSize,
670  sendFaces,
671  sendFacesInfo
672  );
673 
674  // Half2: Determine which faces changed. Use receiveFaces_ ,,
675  nReceiveFaces = getChangedPatchFaces
676  (
677  patch,
678  halfSize,
679  halfSize,
680  receiveFaces,
681  receiveFacesInfo
682  );
683 
684  //Info<< "Half1:" << endl;
685  //writeFaces(nSendFaces, sendFaces, sendFacesInfo, Info);
686  //Info<< endl;
687  //
688  //Info<< "Half2:" << endl;
689  //writeFaces(nReceiveFaces, receiveFaces, receiveFacesInfo, Info);
690  //Info<< endl;
691 
692 
693  // Half1: Adapt wallInfo for leaving domain
694  leaveDomain
695  (
696  patch,
697  nSendFaces,
698  sendFaces,
699  sendFacesInfo
700  );
701  // Half2: Adapt wallInfo for leaving domain
702  leaveDomain
703  (
704  patch,
705  nReceiveFaces,
706  receiveFaces,
707  receiveFacesInfo
708  );
709 
710  // Half1: 'transfer' to other side by offsetting patchFaceI
711  offset(patch, halfSize, nSendFaces, sendFaces);
712 
713  // Half2: 'transfer' to other side
714  offset(patch, -halfSize, nReceiveFaces, receiveFaces);
715 
716  // Apply rotation for non-parallel planes
717  const cyclicPolyPatch& cycPatch =
718  refCast<const cyclicPolyPatch>(patch);
719 
720  if (!cycPatch.parallel())
721  {
722  // sendFaces = received data from half1
723  transform
724  (
725  cycPatch.reverseT(),
726  nSendFaces,
727  sendFacesInfo
728  );
729 
730  // receiveFaces = received data from half2
731  transform
732  (
733  cycPatch.forwardT(),
734  nReceiveFaces,
735  receiveFacesInfo
736  );
737  }
738 
739  if (debug)
740  {
741  Pout<< " Cyclic patch " << patchI << ' ' << patch.name()
742  << " Changed on first half : " << nSendFaces
743  << " Changed on second half : " << nReceiveFaces
744  << endl;
745  }
746 
747  // Half1: Adapt wallInfo for entering domain
748  enterDomain
749  (
750  patch,
751  nSendFaces,
752  sendFaces,
753  sendFacesInfo
754  );
755 
756  // Half2: Adapt wallInfo for entering domain
757  enterDomain
758  (
759  patch,
760  nReceiveFaces,
761  receiveFaces,
762  receiveFacesInfo
763  );
764 
765  // Merge into global storage
766  mergeFaceInfo
767  (
768  patch,
769  nSendFaces,
770  sendFaces,
771  sendFacesInfo,
772  cycPatch.parallel()
773  );
774  // Merge into global storage
775  mergeFaceInfo
776  (
777  patch,
778  nReceiveFaces,
779  receiveFaces,
780  receiveFacesInfo,
781  cycPatch.parallel()
782  );
783 
784  if (debug)
785  {
786  checkCyclic(patch);
787  }
788  }
789  }
790 }
791 
792 
793 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
794 
795 // Set up only. Use setFaceInfo and iterate() to do actual calculation.
796 template <class Type>
798 (
799  const polyMesh& mesh,
800  UList<Type>& allFaceInfo,
801  UList<Type>& allCellInfo
802 )
803 :
804  mesh_(mesh),
805  allFaceInfo_(allFaceInfo),
806  allCellInfo_(allCellInfo),
807  changedFace_(mesh_.nFaces(), false),
808  changedFaces_(mesh_.nFaces()),
809  nChangedFaces_(0),
810  changedCell_(mesh_.nCells(), false),
811  changedCells_(mesh_.nCells()),
812  nChangedCells_(0),
813  hasCyclicPatches_(hasCyclicPatch()),
814  nEvals_(0),
815  nUnvisitedCells_(mesh_.nCells()),
816  nUnvisitedFaces_(mesh_.nFaces()),
817  iter_(0)
818 {}
819 
820 
821 // Iterate, propagating changedFacesInfo across mesh, until no change (or
822 // maxIter reached). Initial cell values specified.
823 template <class Type>
825 (
826  const polyMesh& mesh,
827  const labelList& changedFaces,
828  const List<Type>& changedFacesInfo,
829  UList<Type>& allFaceInfo,
830  UList<Type>& allCellInfo,
831  const label maxIter
832 )
833 :
834  mesh_(mesh),
835  allFaceInfo_(allFaceInfo),
836  allCellInfo_(allCellInfo),
837  changedFace_(mesh_.nFaces(), false),
838  changedFaces_(mesh_.nFaces()),
839  nChangedFaces_(0),
840  changedCell_(mesh_.nCells(), false),
841  changedCells_(mesh_.nCells()),
842  nChangedCells_(0),
843  hasCyclicPatches_(hasCyclicPatch()),
844  nEvals_(0),
845  nUnvisitedCells_(mesh_.nCells()),
846  nUnvisitedFaces_(mesh_.nFaces()),
847  iter_(0)
848 {
849  // Copy initial changed faces data
850  setFaceInfo(changedFaces, changedFacesInfo);
851 
852  // Iterate until nothing changes
853  iterate(maxIter);
854 
855  if ((maxIter > 0) && (iter_ >= maxIter))
856  {
858  (
859  "FaceCellWave<Type>::FaceCellWave"
860  "(const polyMesh&, const labelList&, const List<Type>,"
861  " UList<Type>&, UList<Type>&, const label maxIter)"
862  )
863  << "Maximum number of iterations reached. Increase maxIter." << endl
864  << " maxIter:" << maxIter << endl
865  << " nChangedCells:" << nChangedCells_ << endl
866  << " nChangedFaces:" << nChangedFaces_ << endl
867  << exit(FatalError);
868  }
869 }
870 
871 
872 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
873 
874 
875 template <class Type>
877 {
878  return nUnvisitedCells_;
879 }
880 
881 
882 template <class Type>
884 {
885  return nUnvisitedFaces_;
886 }
887 
888 
889 
890 // Propagate cell to face
891 template <class Type>
893 {
894  const labelList& owner = mesh_.faceOwner();
895  const labelList& neighbour = mesh_.faceNeighbour();
896  label nInternalFaces = mesh_.nInternalFaces();
897 
898  for
899  (
900  label changedFaceI = 0;
901  changedFaceI < nChangedFaces_;
902  changedFaceI++
903  )
904  {
905  label faceI = changedFaces_[changedFaceI];
906  if (!changedFace_[faceI])
907  {
908  FatalErrorIn("FaceCellWave<Type>::faceToCell()")
909  << "Face " << faceI
910  << " not marked as having been changed"
911  << abort(FatalError);
912  }
913 
914 
915  const Type& neighbourWallInfo = allFaceInfo_[faceI];
916 
917  // Evaluate all connected cells
918 
919  // Owner
920  label cellI = owner[faceI];
921  Type& currentWallInfo = allCellInfo_[cellI];
922 
923  if (currentWallInfo != neighbourWallInfo)
924  {
925  updateCell
926  (
927  cellI,
928  faceI,
929  neighbourWallInfo,
930  propagationTol_,
931  currentWallInfo
932  );
933  }
934 
935  // Neighbour. Hack for check if face has neighbour.
936  if (faceI < nInternalFaces)
937  {
938  cellI = neighbour[faceI];
939  Type& currentWallInfo2 = allCellInfo_[cellI];
940 
941  if (currentWallInfo2 != neighbourWallInfo)
942  {
943  updateCell
944  (
945  cellI,
946  faceI,
947  neighbourWallInfo,
948  propagationTol_,
949  currentWallInfo2
950  );
951  }
952  }
953 
954  // Reset status of face
955  changedFace_[faceI] = false;
956  }
957 
958  // Handled all changed faces by now
959  nChangedFaces_ = 0;
960 
961  if (debug)
962  {
963  Pout<< " Changed cells : " << nChangedCells_ << endl;
964  }
965 
966  // Sum nChangedCells over all procs
967  label totNChanged = nChangedCells_;
968 
969  reduce(totNChanged, sumOp<label>());
970 
971  return totNChanged;
972 }
973 
974 
975 // Propagate cell to face
976 template <class Type>
978 {
979  const cellList& cells = mesh_.cells();
980 
981  for
982  (
983  label changedCellI = 0;
984  changedCellI < nChangedCells_;
985  changedCellI++
986  )
987  {
988  label cellI = changedCells_[changedCellI];
989  if (!changedCell_[cellI])
990  {
991  FatalErrorIn("FaceCellWave<Type>::cellToFace()") << "Cell " << cellI
992  << " not marked as having been changed"
993  << abort(FatalError);
994  }
995 
996  const Type& neighbourWallInfo = allCellInfo_[cellI];
997 
998  // Evaluate all connected faces
999 
1000  const labelList& faceLabels = cells[cellI];
1001  forAll(faceLabels, faceLabelI)
1002  {
1003  label faceI = faceLabels[faceLabelI];
1004  Type& currentWallInfo = allFaceInfo_[faceI];
1005 
1006  if (currentWallInfo != neighbourWallInfo)
1007  {
1008  updateFace
1009  (
1010  faceI,
1011  cellI,
1012  neighbourWallInfo,
1013  propagationTol_,
1014  currentWallInfo
1015  );
1016  }
1017  }
1018 
1019  // Reset status of cell
1020  changedCell_[cellI] = false;
1021  }
1022 
1023  // Handled all changed cells by now
1024  nChangedCells_ = 0;
1025 
1026  if (hasCyclicPatches_)
1027  {
1028  // Transfer changed faces across cyclic halves
1029  handleCyclicPatches();
1030  }
1031  if (Pstream::parRun())
1032  {
1033  // Transfer changed faces from neighbouring processors.
1034  handleProcPatches();
1035  }
1036 
1037  if (debug)
1038  {
1039  Pout<< " Changed faces : " << nChangedFaces_ << endl;
1040  }
1041 
1042  // Sum nChangedFaces over all procs
1043  label totNChanged = nChangedFaces_;
1044 
1045  reduce(totNChanged, sumOp<label>());
1046 
1047  return totNChanged;
1048 }
1049 
1050 
1051 // Iterate
1052 template <class Type>
1053 Foam::label Foam::FaceCellWave<Type>::iterate(const label maxIter)
1054 {
1055  if (hasCyclicPatches_)
1056  {
1057  // Transfer changed faces across cyclic halves
1058  handleCyclicPatches();
1059  }
1060  if (Pstream::parRun())
1061  {
1062  // Transfer changed faces from neighbouring processors.
1063  handleProcPatches();
1064  }
1065 
1066  while(iter_ < maxIter)
1067  {
1068  if (debug)
1069  {
1070  Pout<< " Iteration " << iter_ << endl;
1071  }
1072 
1073  nEvals_ = 0;
1074 
1075  label nCells = faceToCell();
1076 
1077  if (debug)
1078  {
1079  Pout<< " Total changed cells : " << nCells << endl;
1080  }
1081 
1082  if (nCells == 0)
1083  {
1084  break;
1085  }
1086 
1087  label nFaces = cellToFace();
1088 
1089  if (debug)
1090  {
1091  Pout<< " Total changed faces : " << nFaces << nl
1092  << " Total evaluations : " << nEvals_ << nl
1093  << " Remaining unvisited cells: " << nUnvisitedCells_ << nl
1094  << " Remaining unvisited faces: " << nUnvisitedFaces_ << endl;
1095  }
1096 
1097  if (nFaces == 0)
1098  {
1099  break;
1100  }
1101 
1102  ++iter_;
1103  }
1104 
1105  return nUnvisitedCells_;
1106 }
1107 
1108 // ************************ vim: set sw=4 sts=4 et: ************************ //