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) 2011-2014 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 "FaceCellWave.H"
27 #include "polyMesh.H"
28 #include "processorPolyPatch.H"
29 #include "cyclicPolyPatch.H"
30 #include "cyclicAMIPolyPatch.H"
31 #include "OPstream.H"
32 #include "IPstream.H"
33 #include "PstreamReduceOps.H"
34 #include "debug.H"
35 #include "typeInfo.H"
36 #include "SubField.H"
37 #include "globalMeshData.H"
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 template<class Type, class TrackingData>
43 
44 template<class Type, class TrackingData>
46 
47 template<class Type, class TrackingData>
49 
50 namespace Foam
51 {
52  //- Combine operator for AMIInterpolation
53  template<class Type, class TrackingData>
54  class combine
55  {
57 
58  const cyclicAMIPolyPatch& patch_;
59 
60  public:
61 
62  combine
63  (
65  const cyclicAMIPolyPatch& patch
66  )
67  :
68  solver_(solver),
69  patch_(patch)
70  {}
71 
72 
73  void operator()
74  (
75  Type& x,
76  const label faceI,
77  const Type& y,
78  const scalar weight
79  ) const
80  {
81  if (y.valid(solver_.data()))
82  {
83  label meshFaceI = -1;
84  if (patch_.owner())
85  {
86  meshFaceI = patch_.start() + faceI;
87  }
88  else
89  {
90  meshFaceI = patch_.neighbPatch().start() + faceI;
91  }
92  x.updateFace
93  (
94  solver_.mesh(),
95  meshFaceI,
96  y,
97  solver_.propagationTol(),
98  solver_.data()
99  );
100  }
101  }
102  };
103 }
104 
105 
106 
107 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
108 
109 // Update info for cellI, at position pt, with information from
110 // neighbouring face/cell.
111 // Updates:
112 // - changedCell_, changedCells_, nChangedCells_,
113 // - statistics: nEvals_, nUnvisitedCells_
114 template<class Type, class TrackingData>
116 (
117  const label cellI,
118  const label neighbourFaceI,
119  const Type& neighbourInfo,
120  const scalar tol,
121  Type& cellInfo
122 )
123 {
124  nEvals_++;
125 
126  bool wasValid = cellInfo.valid(td_);
127 
128  bool propagate =
129  cellInfo.updateCell
130  (
131  mesh_,
132  cellI,
133  neighbourFaceI,
134  neighbourInfo,
135  tol,
136  td_
137  );
138 
139  if (propagate)
140  {
141  if (!changedCell_[cellI])
142  {
143  changedCell_[cellI] = true;
144  changedCells_[nChangedCells_++] = cellI;
145  }
146  }
147 
148  if (!wasValid && cellInfo.valid(td_))
149  {
150  --nUnvisitedCells_;
151  }
152 
153  return propagate;
154 }
155 
156 
157 // Update info for faceI, at position pt, with information from
158 // neighbouring face/cell.
159 // Updates:
160 // - changedFace_, changedFaces_, nChangedFaces_,
161 // - statistics: nEvals_, nUnvisitedFaces_
162 template<class Type, class TrackingData>
164 (
165  const label faceI,
166  const label neighbourCellI,
167  const Type& neighbourInfo,
168  const scalar tol,
169  Type& faceInfo
170 )
171 {
172  nEvals_++;
173 
174  bool wasValid = faceInfo.valid(td_);
175 
176  bool propagate =
177  faceInfo.updateFace
178  (
179  mesh_,
180  faceI,
181  neighbourCellI,
182  neighbourInfo,
183  tol,
184  td_
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(td_))
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, class TrackingData>
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(td_);
222 
223  bool propagate =
224  faceInfo.updateFace
225  (
226  mesh_,
227  faceI,
228  neighbourInfo,
229  tol,
230  td_
231  );
232 
233  if (propagate)
234  {
235  if (!changedFace_[faceI])
236  {
237  changedFace_[faceI] = true;
238  changedFaces_[nChangedFaces_++] = faceI;
239  }
240  }
241 
242  if (!wasValid && faceInfo.valid(td_))
243  {
244  --nUnvisitedFaces_;
245  }
246 
247  return propagate;
248 }
249 
250 
251 // For debugging: check status on both sides of cyclic
252 template<class Type, class TrackingData>
254 (
255  const polyPatch& patch
256 ) const
257 {
258  const cyclicPolyPatch& nbrPatch =
259  refCast<const cyclicPolyPatch>(patch).neighbPatch();
260 
261  forAll(patch, patchFaceI)
262  {
263  label i1 = patch.start() + patchFaceI;
264  label i2 = nbrPatch.start() + patchFaceI;
265 
266  if
267  (
268  !allFaceInfo_[i1].sameGeometry
269  (
270  mesh_,
271  allFaceInfo_[i2],
272  geomTol_,
273  td_
274  )
275  )
276  {
278  (
279  "FaceCellWave<Type, TrackingData>"
280  "::checkCyclic(const polyPatch&)"
281  ) << "problem: i:" << i1 << " otheri:" << i2
282  << " faceInfo:" << allFaceInfo_[i1]
283  << " otherfaceInfo:" << allFaceInfo_[i2]
284  << abort(FatalError);
285  }
286 
287  if (changedFace_[i1] != changedFace_[i2])
288  {
290  (
291  "FaceCellWave<Type, TrackingData>"
292  "::checkCyclic(const polyPatch&)"
293  ) << " problem: i:" << i1 << " otheri:" << i2
294  << " faceInfo:" << allFaceInfo_[i1]
295  << " otherfaceInfo:" << allFaceInfo_[i2]
296  << " changedFace:" << changedFace_[i1]
297  << " otherchangedFace:" << changedFace_[i2]
298  << abort(FatalError);
299  }
300  }
301 }
302 
303 
304 // Check if has cyclic patches
305 template<class Type, class TrackingData>
306 template<class PatchType>
308 {
309  forAll(mesh_.boundaryMesh(), patchI)
310  {
311  if (isA<PatchType>(mesh_.boundaryMesh()[patchI]))
312  {
313  return true;
314  }
315  }
316  return false;
317 }
318 
319 
320 // Copy face information into member data
321 template<class Type, class TrackingData>
323 (
324  const labelList& changedFaces,
325  const List<Type>& changedFacesInfo
326 )
327 {
328  forAll(changedFaces, changedFaceI)
329  {
330  label faceI = changedFaces[changedFaceI];
331 
332  bool wasValid = allFaceInfo_[faceI].valid(td_);
333 
334  // Copy info for faceI
335  allFaceInfo_[faceI] = changedFacesInfo[changedFaceI];
336 
337  // Maintain count of unset faces
338  if (!wasValid && allFaceInfo_[faceI].valid(td_))
339  {
340  --nUnvisitedFaces_;
341  }
342 
343  // Mark faceI as changed, both on list and on face itself.
344 
345  changedFace_[faceI] = true;
346  changedFaces_[nChangedFaces_++] = faceI;
347  }
348 }
349 
350 
351 // Merge face information into member data
352 template<class Type, class TrackingData>
354 (
355  const polyPatch& patch,
356  const label nFaces,
357  const labelList& changedFaces,
358  const List<Type>& changedFacesInfo
359 )
360 {
361  for (label changedFaceI = 0; changedFaceI < nFaces; changedFaceI++)
362  {
363  const Type& neighbourWallInfo = changedFacesInfo[changedFaceI];
364  label patchFaceI = changedFaces[changedFaceI];
365 
366  label meshFaceI = patch.start() + patchFaceI;
367 
368  Type& currentWallInfo = allFaceInfo_[meshFaceI];
369 
370  if (!currentWallInfo.equal(neighbourWallInfo, td_))
371  {
372  updateFace
373  (
374  meshFaceI,
375  neighbourWallInfo,
376  propagationTol_,
377  currentWallInfo
378  );
379  }
380  }
381 }
382 
383 
384 // Construct compact patchFace change arrays for a (slice of a) single patch.
385 // changedPatchFaces in local patch numbering.
386 // Return length of arrays.
387 template<class Type, class TrackingData>
389 (
390  const polyPatch& patch,
391  const label startFaceI,
392  const label nFaces,
393  labelList& changedPatchFaces,
394  List<Type>& changedPatchFacesInfo
395 ) const
396 {
397  label nChangedPatchFaces = 0;
398 
399  for (label i = 0; i < nFaces; i++)
400  {
401  label patchFaceI = i + startFaceI;
402 
403  label meshFaceI = patch.start() + patchFaceI;
404 
405  if (changedFace_[meshFaceI])
406  {
407  changedPatchFaces[nChangedPatchFaces] = patchFaceI;
408  changedPatchFacesInfo[nChangedPatchFaces] = allFaceInfo_[meshFaceI];
409  nChangedPatchFaces++;
410  }
411  }
412  return nChangedPatchFaces;
413 }
414 
415 
416 // Handle leaving domain. Implementation referred to Type
417 template<class Type, class TrackingData>
419 (
420  const polyPatch& patch,
421  const label nFaces,
422  const labelList& faceLabels,
423  List<Type>& faceInfo
424 ) const
425 {
426  const vectorField& fc = mesh_.faceCentres();
427 
428  for (label i = 0; i < nFaces; i++)
429  {
430  label patchFaceI = faceLabels[i];
431 
432  label meshFaceI = patch.start() + patchFaceI;
433  faceInfo[i].leaveDomain(mesh_, patch, patchFaceI, fc[meshFaceI], td_);
434  }
435 }
436 
437 
438 // Handle entering domain. Implementation referred to Type
439 template<class Type, class TrackingData>
441 (
442  const polyPatch& patch,
443  const label nFaces,
444  const labelList& faceLabels,
445  List<Type>& faceInfo
446 ) const
447 {
448  const vectorField& fc = mesh_.faceCentres();
449 
450  for (label i = 0; i < nFaces; i++)
451  {
452  label patchFaceI = faceLabels[i];
453 
454  label meshFaceI = patch.start() + patchFaceI;
455  faceInfo[i].enterDomain(mesh_, patch, patchFaceI, fc[meshFaceI], td_);
456  }
457 }
458 
459 
460 // Transform. Implementation referred to Type
461 template<class Type, class TrackingData>
463 (
464  const tensorField& rotTensor,
465  const label nFaces,
466  List<Type>& faceInfo
467 )
468 {
469  if (rotTensor.size() == 1)
470  {
471  const tensor& T = rotTensor[0];
472 
473  for (label faceI = 0; faceI < nFaces; faceI++)
474  {
475  faceInfo[faceI].transform(mesh_, T, td_);
476  }
477  }
478  else
479  {
480  for (label faceI = 0; faceI < nFaces; faceI++)
481  {
482  faceInfo[faceI].transform(mesh_, rotTensor[faceI], td_);
483  }
484  }
485 }
486 
487 
488 // Offset mesh face. Used for transferring from one cyclic half to the other.
489 template<class Type, class TrackingData>
491 (
492  const polyPatch&,
493  const label cycOffset,
494  const label nFaces,
495  labelList& faces
496 )
497 {
498  for (label faceI = 0; faceI < nFaces; faceI++)
499  {
500  faces[faceI] += cycOffset;
501  }
502 }
503 
504 
505 // Tranfer all the information to/from neighbouring processors
506 template<class Type, class TrackingData>
508 {
509  const globalMeshData& pData = mesh_.globalData();
510 
511  // Which patches are processor patches
512  const labelList& procPatches = pData.processorPatches();
513 
514  // Send all
515 
517 
518  forAll(procPatches, i)
519  {
520  label patchI = procPatches[i];
521 
522  const processorPolyPatch& procPatch =
523  refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
524 
525  // Allocate buffers
526  label nSendFaces;
527  labelList sendFaces(procPatch.size());
528  List<Type> sendFacesInfo(procPatch.size());
529 
530  // Determine which faces changed on current patch
531  nSendFaces = getChangedPatchFaces
532  (
533  procPatch,
534  0,
535  procPatch.size(),
536  sendFaces,
537  sendFacesInfo
538  );
539 
540  // Adapt wallInfo for leaving domain
541  leaveDomain
542  (
543  procPatch,
544  nSendFaces,
545  sendFaces,
546  sendFacesInfo
547  );
548 
549  if (debug & 2)
550  {
551  Pout<< " Processor patch " << patchI << ' ' << procPatch.name()
552  << " communicating with " << procPatch.neighbProcNo()
553  << " Sending:" << nSendFaces
554  << endl;
555  }
556 
557  UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
558  //writeFaces(nSendFaces, sendFaces, sendFacesInfo, toNeighbour);
559  toNeighbour
560  << SubList<label>(sendFaces, nSendFaces)
561  << SubList<Type>(sendFacesInfo, nSendFaces);
562  }
563 
564  pBufs.finishedSends();
565 
566  // Receive all
567 
568  forAll(procPatches, i)
569  {
570  label patchI = procPatches[i];
571 
572  const processorPolyPatch& procPatch =
573  refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
574 
575  // Allocate buffers
576  labelList receiveFaces;
577  List<Type> receiveFacesInfo;
578 
579  {
580  UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
581  fromNeighbour >> receiveFaces >> receiveFacesInfo;
582  }
583 
584  if (debug & 2)
585  {
586  Pout<< " Processor patch " << patchI << ' ' << procPatch.name()
587  << " communicating with " << procPatch.neighbProcNo()
588  << " Receiving:" << receiveFaces.size()
589  << endl;
590  }
591 
592  // Apply transform to received data for non-parallel planes
593  if (!procPatch.parallel())
594  {
595  transform
596  (
597  procPatch.forwardT(),
598  receiveFaces.size(),
599  receiveFacesInfo
600  );
601  }
602 
603  // Adapt wallInfo for entering domain
604  enterDomain
605  (
606  procPatch,
607  receiveFaces.size(),
608  receiveFaces,
609  receiveFacesInfo
610  );
611 
612  // Merge received info
613  mergeFaceInfo
614  (
615  procPatch,
616  receiveFaces.size(),
617  receiveFaces,
618  receiveFacesInfo
619  );
620  }
621 }
622 
623 
624 // Transfer information across cyclic halves.
625 template<class Type, class TrackingData>
627 {
628  forAll(mesh_.boundaryMesh(), patchI)
629  {
630  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
631 
632  if (isA<cyclicPolyPatch>(patch))
633  {
634  const cyclicPolyPatch& nbrPatch =
635  refCast<const cyclicPolyPatch>(patch).neighbPatch();
636 
637  // Allocate buffers
638  label nReceiveFaces;
639  labelList receiveFaces(patch.size());
640  List<Type> receiveFacesInfo(patch.size());
641 
642  // Determine which faces changed
643  nReceiveFaces = getChangedPatchFaces
644  (
645  nbrPatch,
646  0,
647  nbrPatch.size(),
648  receiveFaces,
649  receiveFacesInfo
650  );
651 
652  // Adapt wallInfo for leaving domain
653  leaveDomain
654  (
655  nbrPatch,
656  nReceiveFaces,
657  receiveFaces,
658  receiveFacesInfo
659  );
660 
661  const cyclicPolyPatch& cycPatch =
662  refCast<const cyclicPolyPatch>(patch);
663 
664  if (!cycPatch.parallel())
665  {
666  // received data from other half
667  transform
668  (
669  cycPatch.forwardT(),
670  nReceiveFaces,
671  receiveFacesInfo
672  );
673  }
674 
675  if (debug & 2)
676  {
677  Pout<< " Cyclic patch " << patchI << ' ' << cycPatch.name()
678  << " Changed : " << nReceiveFaces
679  << endl;
680  }
681 
682  // Half2: Adapt wallInfo for entering domain
683  enterDomain
684  (
685  cycPatch,
686  nReceiveFaces,
687  receiveFaces,
688  receiveFacesInfo
689  );
690 
691  // Merge into global storage
692  mergeFaceInfo
693  (
694  cycPatch,
695  nReceiveFaces,
696  receiveFaces,
697  receiveFacesInfo
698  );
699 
700  if (debug)
701  {
702  checkCyclic(cycPatch);
703  }
704  }
705  }
706 }
707 
708 
709 // Transfer information across cyclic halves.
710 template<class Type, class TrackingData>
712 {
713  forAll(mesh_.boundaryMesh(), patchI)
714  {
715  const polyPatch& patch = mesh_.boundaryMesh()[patchI];
716 
717  if (isA<cyclicAMIPolyPatch>(patch))
718  {
719  const cyclicAMIPolyPatch& cycPatch =
720  refCast<const cyclicAMIPolyPatch>(patch);
721 
722  List<Type> receiveInfo;
723 
724  {
725  const cyclicAMIPolyPatch& nbrPatch =
726  refCast<const cyclicAMIPolyPatch>(patch).neighbPatch();
727 
728  // Get nbrPatch data (so not just changed faces)
729  typename List<Type>::subList sendInfo
730  (
731  nbrPatch.patchSlice
732  (
733  allFaceInfo_
734  )
735  );
736 
737  if (!nbrPatch.parallel() || nbrPatch.separated())
738  {
739  // Adapt sendInfo for leaving domain
740  const vectorField::subField fc = nbrPatch.faceCentres();
741  forAll(sendInfo, i)
742  {
743  sendInfo[i].leaveDomain(mesh_, nbrPatch, i, fc[i], td_);
744  }
745  }
746 
747  // Transfer sendInfo to cycPatch
748  combine<Type, TrackingData> cmb(*this, cycPatch);
749 
750  if (cycPatch.applyLowWeightCorrection())
751  {
752  List<Type> defVals
753  (
754  cycPatch.patchInternalList(allCellInfo_)
755  );
756 
757  cycPatch.interpolate(sendInfo, cmb, receiveInfo, defVals);
758  }
759  else
760  {
761  cycPatch.interpolate(sendInfo, cmb, receiveInfo);
762  }
763  }
764 
765  // Apply transform to received data for non-parallel planes
766  if (!cycPatch.parallel())
767  {
768  transform
769  (
770  cycPatch.forwardT(),
771  receiveInfo.size(),
772  receiveInfo
773  );
774  }
775 
776  if (!cycPatch.parallel() || cycPatch.separated())
777  {
778  // Adapt receiveInfo for entering domain
779  const vectorField::subField fc = cycPatch.faceCentres();
780  forAll(receiveInfo, i)
781  {
782  receiveInfo[i].enterDomain(mesh_, cycPatch, i, fc[i], td_);
783  }
784  }
785 
786  // Merge into global storage
787  forAll(receiveInfo, i)
788  {
789  label meshFaceI = cycPatch.start()+i;
790 
791  Type& currentWallInfo = allFaceInfo_[meshFaceI];
792 
793  if
794  (
795  receiveInfo[i].valid(td_)
796  && !currentWallInfo.equal(receiveInfo[i], td_)
797  )
798  {
799  updateFace
800  (
801  meshFaceI,
802  receiveInfo[i],
803  propagationTol_,
804  currentWallInfo
805  );
806  }
807  }
808  }
809  }
810 }
811 
812 
813 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
814 
815 // Set up only. Use setFaceInfo and iterate() to do actual calculation.
816 template<class Type, class TrackingData>
818 (
819  const polyMesh& mesh,
820  UList<Type>& allFaceInfo,
821  UList<Type>& allCellInfo,
822  TrackingData& td
823 )
824 :
825  mesh_(mesh),
826  allFaceInfo_(allFaceInfo),
827  allCellInfo_(allCellInfo),
828  td_(td),
829  changedFace_(mesh_.nFaces(), false),
830  changedFaces_(mesh_.nFaces()),
831  nChangedFaces_(0),
832  changedCell_(mesh_.nCells(), false),
833  changedCells_(mesh_.nCells()),
834  nChangedCells_(0),
835  hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
836  hasCyclicAMIPatches_
837  (
838  returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
839  ),
840  nEvals_(0),
841  nUnvisitedCells_(mesh_.nCells()),
842  nUnvisitedFaces_(mesh_.nFaces())
843 {
844  if
845  (
846  allFaceInfo.size() != mesh_.nFaces()
847  || allCellInfo.size() != mesh_.nCells()
848  )
849  {
851  (
852  "FaceCellWave<Type, TrackingData>::FaceCellWave"
853  "(const polyMesh&, const labelList&, const List<Type>,"
854  " UList<Type>&, UList<Type>&, const label maxIter)"
855  ) << "face and cell storage not the size of mesh faces, cells:"
856  << endl
857  << " allFaceInfo :" << allFaceInfo.size() << endl
858  << " mesh_.nFaces():" << mesh_.nFaces() << endl
859  << " allCellInfo :" << allCellInfo.size() << endl
860  << " mesh_.nCells():" << mesh_.nCells()
861  << exit(FatalError);
862  }
863 }
864 
865 
866 // Iterate, propagating changedFacesInfo across mesh, until no change (or
867 // maxIter reached). Initial cell values specified.
868 template<class Type, class TrackingData>
870 (
871  const polyMesh& mesh,
872  const labelList& changedFaces,
873  const List<Type>& changedFacesInfo,
874  UList<Type>& allFaceInfo,
875  UList<Type>& allCellInfo,
876  const label maxIter,
877  TrackingData& td
878 )
879 :
880  mesh_(mesh),
881  allFaceInfo_(allFaceInfo),
882  allCellInfo_(allCellInfo),
883  td_(td),
884  changedFace_(mesh_.nFaces(), false),
885  changedFaces_(mesh_.nFaces()),
886  nChangedFaces_(0),
887  changedCell_(mesh_.nCells(), false),
888  changedCells_(mesh_.nCells()),
889  nChangedCells_(0),
890  hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
891  hasCyclicAMIPatches_
892  (
893  returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
894  ),
895  nEvals_(0),
896  nUnvisitedCells_(mesh_.nCells()),
897  nUnvisitedFaces_(mesh_.nFaces())
898 {
899  if
900  (
901  allFaceInfo.size() != mesh_.nFaces()
902  || allCellInfo.size() != mesh_.nCells()
903  )
904  {
906  (
907  "FaceCellWave<Type, TrackingData>::FaceCellWave"
908  "(const polyMesh&, const labelList&, const List<Type>,"
909  " UList<Type>&, UList<Type>&, const label maxIter)"
910  ) << "face and cell storage not the size of mesh faces, cells:"
911  << endl
912  << " allFaceInfo :" << allFaceInfo.size() << endl
913  << " mesh_.nFaces():" << mesh_.nFaces() << endl
914  << " allCellInfo :" << allCellInfo.size() << endl
915  << " mesh_.nCells():" << mesh_.nCells()
916  << exit(FatalError);
917  }
918 
919  // Copy initial changed faces data
920  setFaceInfo(changedFaces, changedFacesInfo);
921 
922  // Iterate until nothing changes
923  label iter = iterate(maxIter);
924 
925  if ((maxIter > 0) && (iter >= maxIter))
926  {
928  (
929  "FaceCellWave<Type, TrackingData>::FaceCellWave"
930  "(const polyMesh&, const labelList&, const List<Type>,"
931  " UList<Type>&, UList<Type>&, const label maxIter)"
932  )
933  << "Maximum number of iterations reached. Increase maxIter." << endl
934  << " maxIter:" << maxIter << endl
935  << " nChangedCells:" << nChangedCells_ << endl
936  << " nChangedFaces:" << nChangedFaces_ << endl
937  << exit(FatalError);
938  }
939 }
940 
941 
942 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
943 
944 
945 template<class Type, class TrackingData>
947 {
948  return nUnvisitedCells_;
949 }
950 
951 
952 template<class Type, class TrackingData>
954 {
955  return nUnvisitedFaces_;
956 }
957 
958 
959 
960 // Propagate cell to face
961 template<class Type, class TrackingData>
963 {
964  const labelList& owner = mesh_.faceOwner();
965  const labelList& neighbour = mesh_.faceNeighbour();
966  label nInternalFaces = mesh_.nInternalFaces();
967 
968  for
969  (
970  label changedFaceI = 0;
971  changedFaceI < nChangedFaces_;
972  changedFaceI++
973  )
974  {
975  label faceI = changedFaces_[changedFaceI];
976  if (!changedFace_[faceI])
977  {
978  FatalErrorIn("FaceCellWave<Type, TrackingData>::faceToCell()")
979  << "Face " << faceI
980  << " not marked as having been changed"
981  << abort(FatalError);
982  }
983 
984 
985  const Type& neighbourWallInfo = allFaceInfo_[faceI];
986 
987  // Evaluate all connected cells
988 
989  // Owner
990  label cellI = owner[faceI];
991  Type& currentWallInfo = allCellInfo_[cellI];
992 
993  if (!currentWallInfo.equal(neighbourWallInfo, td_))
994  {
995  updateCell
996  (
997  cellI,
998  faceI,
999  neighbourWallInfo,
1000  propagationTol_,
1001  currentWallInfo
1002  );
1003  }
1004 
1005  // Neighbour.
1006  if (faceI < nInternalFaces)
1007  {
1008  cellI = neighbour[faceI];
1009  Type& currentWallInfo2 = allCellInfo_[cellI];
1010 
1011  if (!currentWallInfo2.equal(neighbourWallInfo, td_))
1012  {
1013  updateCell
1014  (
1015  cellI,
1016  faceI,
1017  neighbourWallInfo,
1018  propagationTol_,
1019  currentWallInfo2
1020  );
1021  }
1022  }
1023 
1024  // Reset status of face
1025  changedFace_[faceI] = false;
1026  }
1027 
1028  // Handled all changed faces by now
1029  nChangedFaces_ = 0;
1030 
1031  if (debug & 2)
1032  {
1033  Pout<< " Changed cells : " << nChangedCells_ << endl;
1034  }
1035 
1036  // Sum nChangedCells over all procs
1037  label totNChanged = nChangedCells_;
1038 
1039  reduce(totNChanged, sumOp<label>());
1040 
1041  return totNChanged;
1042 }
1043 
1044 
1045 // Propagate cell to face
1046 template<class Type, class TrackingData>
1048 {
1049  const cellList& cells = mesh_.cells();
1050 
1051  for
1052  (
1053  label changedCellI = 0;
1054  changedCellI < nChangedCells_;
1055  changedCellI++
1056  )
1057  {
1058  label cellI = changedCells_[changedCellI];
1059  if (!changedCell_[cellI])
1060  {
1061  FatalErrorIn("FaceCellWave<Type, TrackingData>::cellToFace()")
1062  << "Cell " << cellI << " not marked as having been changed"
1063  << abort(FatalError);
1064  }
1065 
1066  const Type& neighbourWallInfo = allCellInfo_[cellI];
1067 
1068  // Evaluate all connected faces
1069 
1070  const labelList& faceLabels = cells[cellI];
1071  forAll(faceLabels, faceLabelI)
1072  {
1073  label faceI = faceLabels[faceLabelI];
1074  Type& currentWallInfo = allFaceInfo_[faceI];
1075 
1076  if (!currentWallInfo.equal(neighbourWallInfo, td_))
1077  {
1078  updateFace
1079  (
1080  faceI,
1081  cellI,
1082  neighbourWallInfo,
1083  propagationTol_,
1084  currentWallInfo
1085  );
1086  }
1087  }
1088 
1089  // Reset status of cell
1090  changedCell_[cellI] = false;
1091  }
1092 
1093  // Handled all changed cells by now
1094  nChangedCells_ = 0;
1095 
1096  if (hasCyclicPatches_)
1097  {
1098  // Transfer changed faces across cyclic halves
1099  handleCyclicPatches();
1100  }
1101 
1102  if (hasCyclicAMIPatches_)
1103  {
1104  handleAMICyclicPatches();
1105  }
1106 
1107  if (Pstream::parRun())
1108  {
1109  // Transfer changed faces from neighbouring processors.
1110  handleProcPatches();
1111  }
1112 
1113  if (debug & 2)
1114  {
1115  Pout<< " Changed faces : " << nChangedFaces_ << endl;
1116  }
1117 
1118  // Sum nChangedFaces over all procs
1119  label totNChanged = nChangedFaces_;
1120 
1121  reduce(totNChanged, sumOp<label>());
1122 
1123  return totNChanged;
1124 }
1125 
1126 
1127 // Iterate
1128 template<class Type, class TrackingData>
1130 {
1131  if (hasCyclicPatches_)
1132  {
1133  // Transfer changed faces across cyclic halves
1134  handleCyclicPatches();
1135  }
1136 
1137  if (hasCyclicAMIPatches_)
1138  {
1139  handleAMICyclicPatches();
1140  }
1141 
1142  if (Pstream::parRun())
1143  {
1144  // Transfer changed faces from neighbouring processors.
1145  handleProcPatches();
1146  }
1147 
1148  label iter = 0;
1149 
1150  while (iter < maxIter)
1151  {
1152  if (debug)
1153  {
1154  Info<< " Iteration " << iter << endl;
1155  }
1156 
1157  nEvals_ = 0;
1158 
1159  label nCells = faceToCell();
1160 
1161  if (debug)
1162  {
1163  Info<< " Total changed cells : " << nCells << endl;
1164  }
1165 
1166  if (nCells == 0)
1167  {
1168  break;
1169  }
1170 
1171  label nFaces = cellToFace();
1172 
1173  if (debug)
1174  {
1175  Info<< " Total changed faces : " << nFaces << nl
1176  << " Total evaluations : " << nEvals_ << nl
1177  << " Remaining unvisited cells: " << nUnvisitedCells_ << nl
1178  << " Remaining unvisited faces: " << nUnvisitedFaces_ << endl;
1179  }
1180 
1181  if (nFaces == 0)
1182  {
1183  break;
1184  }
1185 
1186  ++iter;
1187  }
1188 
1189  return iter;
1190 }
1191 
1192 
1193 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual bool owner() const
Does this side own the patch?
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:380
combine(FaceCellWave< Type, TrackingData > &solver, const cyclicAMIPolyPatch &patch)
Definition: FaceCellWave.C:63
virtual const tensorField & forwardT() const
Return face transformation tensor.
Pre-declare related SubField type.
Definition: Field.H:61
Buffers for inter-processor communications streams (UOPstream, UIPstream).
Cyclic patch for Arbitrary Mesh Interface (AMI)
virtual bool separated() const
Are the planes separated.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void setFaceInfo(const labelList &changedFaces, const List< Type > &changedFacesInfo)
Set initial changed faces.
Definition: FaceCellWave.C:323
const word & name() const
Return name.
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)
A topoSetSource to select cells based on usage in faces.
Definition: faceToCell.H:49
label getUnsetCells() const
Get number of unvisited cells, i.e. cells that were not (yet)
Definition: FaceCellWave.C:946
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
messageStream Info
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:278
Namespace for OpenFOAM.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:284
A topoSetSource to select a faceSet from cells.
Definition: cellToFace.H:53
virtual const cyclicAMIPolyPatch & neighbPatch() const
Return a reference to the neighbour patch.
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Combine operator for AMIInterpolation.
Definition: FaceCellWave.C:54
Neighbour processor patch.
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition: FaceCellWave.H:74
virtual bool parallel() const
Are the cyclic planes parallel.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
#define forAll(list, i)
Definition: UList.H:421
label cellToFace()
Propagate from cell to face. Returns total number of faces.
tmp< Field< Type > > interpolate(const Field< Type > &fld, const UList< Type > &defaultValues=UList< Type >()) const
Interpolate field.
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
Cyclic plane patch.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
const List< T >::subList patchSlice(const UList< T > &l) const
Slice list to patch.
Definition: polyPatch.H:333
bool applyLowWeightCorrection() const
Return true if applying the low weight correction.
const cellShapeList & cells
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
Definition: cellInfo.H:54
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label faceToCell()
Propagate from face to cell. Returns total number of cells.
Definition: FaceCellWave.C:962
error FatalError
scalar y
const polyMesh & mesh() const
Access mesh.
Definition: FaceCellWave.H:312
label getUnsetFaces() const
Get number of unvisited faces.
Definition: FaceCellWave.C:953
int neighbProcNo() const
Return neigbour processor number.
const TrackingData & data() const
Additional data to be passed into container.
Definition: FaceCellWave.H:306
A List obtained as a section of another List.
Definition: SubList.H:53
static scalar propagationTol()
Access to tolerance.
Definition: FaceCellWave.H:250
label iterate(const label maxIter)
Iterate until no changes or maxIter reached. Returns actual.
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
const UIndirectList< T > patchInternalList(const UList< T > &internalValues) const
Extract face cell data.
Definition: polyPatch.H:324
const labelList & processorPatches() const
Return list of processor patch labels.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53