FreeFOAM The Cross-Platform CFD Toolkit
referredCellList.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) 2008-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 "referredCellList.H"
30 #include <OpenFOAM/Time.H>
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 void Foam::referredCellList::buildReferredCellList
38 (
39  bool pointPointListBuild
40 )
41 {
42  Info << nl << "Building list of referred interaction neighbours" << endl;
43 
44  const polyMesh& mesh(il_.mesh());
45 
46  DynamicList<referredCell> referredInteractionList;
47 
48  // realCellsWithinRCutMaxOfAnyReferringPatch
49  DynamicList<label> rCellsWRRP;
50 
51  // realFacesWithinRCutMaxOfAnyReferringPatch
52  DynamicList<label> rFacesWRRP;
53 
54  // realEdgesWithinRCutMaxOfAnyReferringPatch
55  DynamicList<label> rEdgesWRRP;
56 
57  // realPointsWithinRCutMaxOfAnyReferringPatch
58  DynamicList<label> rPointsWRRP;
59 
60  labelListList processorPatchSegmentMapping
61  (
62  mesh.globalData().processorPatches().size()
63  );
64 
65  List<vectorList> allNeighbourFaceCentres
66  (
67  mesh.globalData().processorPatches().size()
68  );
69 
70  List<vectorList> allNeighbourFaceAreas
71  (
72  mesh.globalData().processorPatches().size()
73  );
74 
75  label nUndecomposedPatches = 0;
76 
77  if (Pstream::parRun())
78  {
79  dictionary patchDictionary;
80 
81  DynamicList<word> patchNames;
82 
83  Time undecomposedTime
84  (
86  mesh.time().rootPath(),
87  mesh.time().caseName().path()
88  );
89 
90  IOobject undecomposedBoundaryHeader
91  (
92  "boundary",
93  undecomposedTime.constant(),
95  undecomposedTime,
98  false
99  );
100 
101  if (undecomposedBoundaryHeader.headerOk())
102  {
103  polyBoundaryMeshEntries undecomposedPatchEntries
104  (
105  undecomposedBoundaryHeader
106  );
107 
108  forAll(undecomposedPatchEntries, patchi)
109  {
110  patchNames.append
111  (
112  undecomposedPatchEntries[patchi].keyword()
113  );
114 
115  patchDictionary.add
116  (
117  undecomposedPatchEntries[patchi]
118  );
119  }
120  }
121  else
122  {
123  FatalErrorIn ("referredCellList.C")
124  << nl << "unable to read undecomposed boundary file from "
125  << "constant/polyMesh" << nl
126  << abort(FatalError);
127  }
128 
129  labelIOList faceProcAddressing
130  (
131  IOobject
132  (
133  "faceProcAddressing",
134  mesh.time().constant(),
136  mesh,
137  IOobject::MUST_READ,
138  IOobject::NO_WRITE,
139  false
140  )
141  );
142 
143  labelList procPatches(mesh.globalData().processorPatches());
144 
145  nUndecomposedPatches = patchNames.size();
146 
147  // processorPatchSegmentMapping works by mapping the original patch and
148  // half that a face on a processor patch was on before decomposition.
149  // This creates a patch segment for each half of each original (cyclic)
150  // patch which can be assessed separately. There are n =
151  // patchNames.size() original patches, k = 0 to n-1. The mapping is:
152  // value = 0: originally an internal face value = k, was originally on
153  // the on the patch k-1, 1st half value = -k, was originally on the on
154  // the patch k-1, 2nd half
155 
156  forAll(procPatches,pP)
157  {
158  const processorPolyPatch& patch = refCast<const processorPolyPatch>
159  (
160  mesh.boundaryMesh()[procPatches[pP]]
161  );
162 
163  labelList& procPatchSegMap = processorPatchSegmentMapping[pP];
164 
165  procPatchSegMap.setSize(patch.size());
166 
167  forAll (patch, pI)
168  {
169  label decomposedMeshFace = patch.start() + pI;
170 
171  label faceProcAdd = faceProcAddressing[decomposedMeshFace];
172 
173  label globalFace = abs(faceProcAdd)-1;
174 
175  label minStart = -1;
176 
177  forAll(patchNames, pN)
178  {
179  if (patchDictionary.found(patchNames[pN]))
180  {
181  const dictionary& patchDict =
182  patchDictionary.subDict(patchNames[pN]);
183 
184  label startFace
185  (
186  readLabel
187  (
188  patchDict.lookup("startFace")
189  )
190  );
191 
192  label nFaces(readLabel(patchDict.lookup("nFaces")));
193 
194  if (minStart < 0 || startFace < minStart)
195  {
196  minStart = startFace;
197  }
198 
199  if
200  (
201  globalFace >= startFace
202  && globalFace < startFace + nFaces/2
203  )
204  {
205  procPatchSegMap[pI] = pN + 1;
206  }
207  else if
208  (
209  globalFace >= startFace + nFaces/2
210  && globalFace < startFace + nFaces
211  )
212  {
213  procPatchSegMap[pI] = -(pN + 1);
214  }
215  }
216  }
217 
218  if (globalFace < minStart)
219  {
220  procPatchSegMap[pI] = 0;
221  }
222  }
223  }
224 
225  forAll(procPatches,pP)
226  {
227  const processorPolyPatch& patch = refCast<const processorPolyPatch>
228  (
229  mesh.boundaryMesh()[procPatches[pP]]
230  );
231 
232  {
233  OPstream toNeighbProc
234  (
236  patch.neighbProcNo()
237  );
238 
239  toNeighbProc << patch.faceCentres() << patch.faceAreas();
240  }
241  }
242 
243  forAll(procPatches,pP)
244  {
245  const processorPolyPatch& patch = refCast<const processorPolyPatch>
246  (
247  mesh.boundaryMesh()[procPatches[pP]]
248  );
249 
250  vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
251 
252  neighbFaceCentres.setSize(patch.size());
253 
254  vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
255 
256  neighbFaceAreas.setSize(patch.size());
257 
258  {
259  IPstream fromNeighbProc
260  (
262  patch.neighbProcNo()
263  );
264 
265  fromNeighbProc >> neighbFaceCentres >> neighbFaceAreas;
266  }
267  }
268 
269  // *************************************************************
270  // Tests that all processor patch segments from different
271  // original patches prior to decomposition produce the same
272  // transform. Check before 1st iteration.
273  // *************************************************************
274 
275  forAll(procPatches,pP)
276  {
277  const processorPolyPatch& patch = refCast<const processorPolyPatch>
278  (
279  mesh.boundaryMesh()[procPatches[pP]]
280  );
281 
282  const vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
283 
284  const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
285 
286  label nUP;
287 
288  for
289  (
290  nUP = -nUndecomposedPatches;
291  nUP <= nUndecomposedPatches;
292  nUP++
293  )
294  {
295  DynamicList<vector> refOff;
296 
297  DynamicList<tensor> refTrans;
298 
299  forAll (patch, faceI)
300  {
301  if (processorPatchSegmentMapping[pP][faceI] == nUP)
302  {
303  referredCell testRefCell
304  (
305  mesh,
306  -1,
307  -1,
308  patch.faceCentres()[faceI],
309  neighbFaceCentres[faceI],
310  patch.faceNormals()[faceI],
311  neighbFaceAreas[faceI]
312  /(mag(neighbFaceAreas[faceI]) + VSMALL)
313  );
314 
315  refOff.append(testRefCell.offset());
316 
317  refTrans.append(testRefCell.rotation());
318  }
319  }
320 
321  refOff.shrink();
322 
323  refTrans.shrink();
324 
325  if (refOff.size())
326  {
327  if
328  (
329  sum(mag(refOff-refOff[0]))/refOff.size()
331  || sum(mag(refTrans-refTrans[0]))/refTrans.size()
333  )
334  {
335  FatalErrorIn ("referredCellList.C")
336  << nl << "Face pairs on patch "
337  << patch.name()
338  << ", segment " << patchNames[nUP]
339  << " do not give the same referring "
340  << " transformations to within tolerance of "
342  << " Referring offsets:" << refOff << nl
343  << " Average sum of mag difference: "
344  << sum(mag(refOff-refOff[0]))/refOff.size() << nl
345  << " Referring transforms:" << refTrans << nl
346  << " Average sum of mag difference: "
347  << sum(mag(refTrans-refTrans[0]))/refTrans.size()
348  << nl << abort(FatalError);
349  }
350  }
351  }
352  }
353  }
354 
355  label cellsReferredThisIteration = 1;
356 
357  label iterationNo = 0;
358 
359  while (cellsReferredThisIteration)
360  {
361  label refIntListStartSize = referredInteractionList.size();
362 
363  forAll (mesh.boundaryMesh(), patchI)
364  {
365  // Treat local cyclics on each processor before processor
366  // boundaries. Separate treatment allows the serial version to run
367  // transparently.
368 
369  if (isA<cyclicPolyPatch>(mesh.boundaryMesh()[patchI]))
370  {
371  const cyclicPolyPatch& patch = refCast<const cyclicPolyPatch>
372  (
373  mesh.boundaryMesh()[patchI]
374  );
375 
376  if (patch.size())
377  {
378  if (iterationNo == 0)
379  {
380  // Tests that all combinations of face pairs produce the
381  // same transform. Only check on 1st iteration.
382 
383  label faceL;
384  // A face in the 1st half of the patch
385 
386  label faceM;
387  // Face corresponding to faceL in the 2nd half of the
388  // patch. Assumes correct face ordering.
389 
390  vectorList refOff(patch.size()/2);
391 
392  List<tensor> refTrans(patch.size()/2);
393 
394  for
395  (
396  faceL = 0, faceM = patch.size()/2;
397  faceL < patch.size()/2;
398  faceL++, faceM++
399  )
400  {
401  referredCell testRefCell
402  (
403  mesh,
404  -1,
405  -1,
406  patch.faceCentres()[faceL],
407  patch.faceCentres()[faceM],
408  patch.faceNormals()[faceL],
409  patch.faceNormals()[faceM]
410  );
411 
412  refOff[faceL] = testRefCell.offset();
413 
414  refTrans[faceL] = testRefCell.rotation();
415  }
416 
417  if
418  (
419  sum(mag(refOff - refOff[0]))/(patch.size()/2)
421  || sum(mag(refTrans - refTrans[0]))/(patch.size()/2)
423  )
424  {
425  FatalErrorIn ("referredCellList.C")
426  << nl << "Face pairs on patch "
427  << patch.name()
428  << " do not give the same referring "
429  << " transformations to within tolerance of "
431  << " Referring offsets:" << refOff << nl
432  << " Average sum of mag difference: "
433  << sum(mag(refOff-refOff[0]))/refOff.size()
434  << nl
435  << " Referring transforms:" << refTrans << nl
436  << " Average sum of mag difference: "
437  << sum(mag(refTrans-refTrans[0]))
438  /refTrans.size()
439  << nl << abort(FatalError);
440  }
441  }
442 
443  // *********************************************************
444  // 1st half of face list - 1st side of boundary
445  // *********************************************************
446 
447  label faceI;
448 
449  DynamicList<label> meshFacesOnThisSegment;
450 
451  for (faceI = 0; faceI < patch.size()/2; faceI++)
452  {
453  // unable to use the normal accessors for the polyPatch
454  // because points on separate halves need used
455  // separately
456 
457  meshFacesOnThisSegment.append(faceI + patch.start());
458  }
459 
460  meshFacesOnThisSegment.shrink();
461 
462  DynamicList<label> meshEdgesOnThisSegment;
463 
464  DynamicList<label> meshPointsOnThisSegment;
465 
466  forAll(meshFacesOnThisSegment, mFOTS)
467  {
468  const label segFace = meshFacesOnThisSegment[mFOTS];
469 
470  const labelList& faceEdges = mesh.faceEdges()[segFace];
471 
472  forAll (faceEdges, fE)
473  {
474  const label faceEdge(faceEdges[fE]);
475 
476  if
477  (
478  findIndex
479  (
480  meshEdgesOnThisSegment,
481  faceEdge
482  ) == -1
483  )
484  {
485  meshEdgesOnThisSegment.append(faceEdge);
486  }
487  }
488 
489  const face& facePoints(mesh.faces()[segFace]);
490 
491  forAll (facePoints, fP)
492  {
493  const label facePoint(facePoints[fP]);
494 
495  if
496  (
497  findIndex
498  (
499  meshPointsOnThisSegment,
500  facePoint
501  )
502  ==
503  -1
504  )
505  {
506  meshPointsOnThisSegment.append(facePoint);
507  }
508  }
509  }
510 
511  meshEdgesOnThisSegment.shrink();
512 
513  meshPointsOnThisSegment.shrink();
514 
515  if (iterationNo == 0)
516  {
517  // Assessing real cells in range is only required on
518  // the 1st iteration because they do not change from
519  // iteration to iteration.
520 
521  labelList realCellsFoundInRange
522  (
524  (
525  meshFacesOnThisSegment,
526  meshEdgesOnThisSegment,
527  meshPointsOnThisSegment
528  )
529  );
530 
531  forAll(realCellsFoundInRange,cFIR)
532  {
533  const label realCell = realCellsFoundInRange[cFIR];
534 
535  referredCell cellToRefer
536  (
537  mesh,
539  realCell,
540  patch.faceCentres()[0],
541  patch.faceCentres()[patch.size()/2],
542  patch.faceNormals()[0],
543  patch.faceNormals()[patch.size()/2]
544  );
545 
546  // Test all existing referred and real cells to
547  // check duplicates are not being made or cells
548  // aren't being referred back onto themselves
549 
550  bool addCellToRefer = true;
551 
552  // Check if cellToRefer is an existing referred cell
553 
554  forAll(referredInteractionList, rIL)
555  {
556  if
557  (
558  cellToRefer.duplicate
559  (
560  referredInteractionList[rIL]
561  )
562  )
563  {
564  addCellToRefer = false;
565 
566  break;
567  }
568  }
569 
570  // Check for cellToRefer being referred back
571  // ontop of a real cell
572 
573  if
574  (
575  cellToRefer.duplicate
576  (
578  mesh.nCells()
579  )
580  )
581  {
582  addCellToRefer = false;
583  }
584 
585  if (addCellToRefer)
586  {
587  referredInteractionList.append(cellToRefer);
588  }
589 
590  // add real cells found in range of cyclic patch
591  // to whole mesh list
592 
593  if (findIndex (rCellsWRRP, realCell) == -1)
594  {
595  rCellsWRRP.append(realCell);
596  }
597  }
598  }
599 
600  referredInteractionList.shrink();
601 
602  labelList referredCellsFoundInRange
603  (
605  (
606  referredInteractionList,
607  meshFacesOnThisSegment,
608  meshEdgesOnThisSegment,
609  meshPointsOnThisSegment
610  )
611  );
612 
613  forAll(referredCellsFoundInRange,cFIR)
614  {
615  referredCell& existingRefCell =
616  referredInteractionList
617  [
618  referredCellsFoundInRange[cFIR]
619  ];
620 
621  referredCell cellToReRefer =
622  existingRefCell.reRefer
623  (
624  patch.faceCentres()[0],
625  patch.faceCentres()[patch.size()/2],
626  patch.faceNormals()[0],
627  patch.faceNormals()[patch.size()/2]
628  );
629 
630  // Test all existing referred and real cells to check
631  // duplicates are not being made or cells aren't being
632  // referred back onto themselves
633 
634  bool addCellToReRefer = true;
635 
636  // Check if cellToRefer is an existing referred cell
637 
638  forAll(referredInteractionList, rIL)
639  {
640  if
641  (
642  cellToReRefer.duplicate
643  (
644  referredInteractionList[rIL]
645  )
646  )
647  {
648  addCellToReRefer = false;
649 
650  break;
651  }
652  }
653 
654  // Check for cellToRefer being referred back
655  // ontop of a real cell
656 
657  if
658  (
659  cellToReRefer.duplicate
660  (
662  mesh.nCells()
663  )
664  )
665  {
666  addCellToReRefer = false;
667  }
668 
669  if (addCellToReRefer)
670  {
671  referredInteractionList.append(cellToReRefer);
672  }
673  }
674 
675  // *********************************************************
676  // 2nd half of face list - 2nd side of boundary
677  // *********************************************************
678 
679  meshFacesOnThisSegment.clear();
680 
681  for (faceI = patch.size()/2; faceI < patch.size(); faceI++)
682  {
683  // unable to use the normal accessors for the polyPatch
684  // because points on separate halves need used
685  // separately
686 
687  meshFacesOnThisSegment.append(faceI + patch.start());
688  }
689 
690  meshFacesOnThisSegment.shrink();
691 
692  meshEdgesOnThisSegment.clear();
693 
694  meshPointsOnThisSegment.clear();
695 
696  forAll(meshFacesOnThisSegment, mFOTS)
697  {
698  const label segFace = meshFacesOnThisSegment[mFOTS];
699 
700  const labelList& faceEdges = mesh.faceEdges()[segFace];
701 
702  forAll (faceEdges, fE)
703  {
704  const label faceEdge(faceEdges[fE]);
705 
706  if
707  (
708  findIndex
709  (
710  meshEdgesOnThisSegment,
711  faceEdge
712  )
713  ==
714  -1
715  )
716  {
717  meshEdgesOnThisSegment.append(faceEdge);
718  }
719  }
720 
721  const face& facePoints(mesh.faces()[segFace]);
722 
723  forAll (facePoints, fP)
724  {
725  const label facePoint(facePoints[fP]);
726 
727  if
728  (
729  findIndex
730  (
731  meshPointsOnThisSegment,
732  facePoint
733  )
734  ==
735  -1
736  )
737  {
738  meshPointsOnThisSegment.append(facePoint);
739  }
740  }
741  }
742 
743  meshEdgesOnThisSegment.shrink();
744 
745  meshPointsOnThisSegment.shrink();
746 
747  if (iterationNo == 0)
748  {
749  // Assessing real cells in range is only required on
750  // the 1st iteration because they do not change from
751  // iteration to iteration.
752 
753  labelList realCellsFoundInRange
754  (
756  (
757  meshFacesOnThisSegment,
758  meshEdgesOnThisSegment,
759  meshPointsOnThisSegment
760  )
761  );
762 
763  forAll(realCellsFoundInRange,cFIR)
764  {
765  const label realCell = realCellsFoundInRange[cFIR];
766 
767  referredCell cellToRefer
768  (
769  mesh,
771  realCell,
772  patch.faceCentres()[patch.size()/2],
773  patch.faceCentres()[0],
774  patch.faceNormals()[patch.size()/2],
775  patch.faceNormals()[0]
776  );
777 
778  // Test all existing referred and real cells to
779  // check duplicates are not being made or cells
780  // aren't being referred back onto themselves
781 
782  bool addCellToRefer = true;
783 
784  // Check if cellToRefer is an existing referred cell
785 
786  forAll(referredInteractionList, rIL)
787  {
788  if
789  (
790  cellToRefer.duplicate
791  (
792  referredInteractionList[rIL]
793  )
794  )
795  {
796  addCellToRefer = false;
797 
798  break;
799  }
800  }
801 
802  // Check for cellToRefer being referred back
803  // ontop of a real cell
804 
805  if
806  (
807  cellToRefer.duplicate
808  (
810  mesh.nCells()
811  )
812  )
813  {
814  addCellToRefer = false;
815  }
816 
817  if (addCellToRefer)
818  {
819  referredInteractionList.append(cellToRefer);
820  }
821 
822  // add real cells found in range of cyclic patch
823  // to whole mesh list
824 
825  if (findIndex (rCellsWRRP, realCell) == -1)
826  {
827  rCellsWRRP.append(realCell);
828  }
829  }
830  }
831 
832  referredInteractionList.shrink();
833 
834  referredCellsFoundInRange =
836  (
837  referredInteractionList,
838  meshFacesOnThisSegment,
839  meshEdgesOnThisSegment,
840  meshPointsOnThisSegment
841  );
842 
843  forAll(referredCellsFoundInRange,cFIR)
844  {
845  referredCell& existingRefCell =
846  referredInteractionList
847  [
848  referredCellsFoundInRange[cFIR]
849  ];
850 
851  referredCell cellToReRefer =
852  existingRefCell.reRefer
853  (
854  patch.faceCentres()[patch.size()/2],
855  patch.faceCentres()[0],
856  patch.faceNormals()[patch.size()/2],
857  patch.faceNormals()[0]
858  );
859 
860  // Test all existing referred and real cells to check
861  // duplicates are not being made or cells aren't being
862  // referred back onto themselves
863 
864  bool addCellToReRefer = true;
865 
866  // Check if cellToRefer is an existing referred cell
867 
868  forAll(referredInteractionList, rIL)
869  {
870  if
871  (
872  cellToReRefer.duplicate
873  (
874  referredInteractionList[rIL]
875  )
876  )
877  {
878  addCellToReRefer = false;
879 
880  break;
881  }
882  }
883 
884  // Check for cellToRefer being referred back
885  // ontop of a real cell
886 
887  if
888  (
889  cellToReRefer.duplicate
890  (
892  mesh.nCells()
893  )
894  )
895  {
896  addCellToReRefer = false;
897  }
898 
899  if (addCellToReRefer)
900  {
901  referredInteractionList.append(cellToReRefer);
902  }
903  }
904  }
905  }
906  }
907 
908  if (Pstream::parRun())
909  {
910  labelList procPatches(mesh.globalData().processorPatches());
911 
912  forAll(procPatches,pP)
913  {
914  const processorPolyPatch& patch =
915  refCast<const processorPolyPatch>
916  (
917  mesh.boundaryMesh()[procPatches[pP]]
918  );
919 
920  DynamicList<referredCell> referredCellsToTransfer;
921 
922  const vectorList& neighbFaceCentres =
923  allNeighbourFaceCentres[pP];
924 
925  const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
926 
927  label nUP;
928 
929  for
930  (
931  nUP = -nUndecomposedPatches;
932  nUP <= nUndecomposedPatches;
933  nUP++
934  )
935  {
936  // faceT is used to specify one face on this patch segment
937  // that will be used to calculate the transformation values.
938  // All faces are guaranteed to produce the same transform
939  // because of the checks carried out at the start of the
940  // function. Setting to -1 until the 1st face on this
941  // segment is found.
942 
943  label faceT = -1;
944 
945  DynamicList<label> meshFacesOnThisSegment;
946 
947  forAll (patch, faceI)
948  {
949  if (processorPatchSegmentMapping[pP][faceI] == nUP)
950  {
951  if (faceT == -1)
952  {
953  faceT = faceI;
954  }
955 
956  meshFacesOnThisSegment.append
957  (
958  faceI + patch.start()
959  );
960  }
961  }
962 
963  meshFacesOnThisSegment.shrink();
964 
965  DynamicList<label> meshEdgesOnThisSegment;
966 
967  DynamicList<label> meshPointsOnThisSegment;
968 
969  forAll(meshFacesOnThisSegment, mFOTS)
970  {
971  const label segFace = meshFacesOnThisSegment[mFOTS];
972 
973  const labelList& faceEdges = mesh.faceEdges()[segFace];
974 
975  forAll (faceEdges, fE)
976  {
977  const label faceEdge(faceEdges[fE]);
978 
979  if
980  (
981  findIndex
982  (
983  meshEdgesOnThisSegment,
984  faceEdge
985  )
986  ==
987  -1
988  )
989  {
990  meshEdgesOnThisSegment.append(faceEdge);
991  }
992  }
993 
994  const face& facePoints(mesh.faces()[segFace]);
995 
996  forAll (facePoints, fP)
997  {
998  const label facePoint(facePoints[fP]);
999 
1000  if
1001  (
1002  findIndex
1003  (
1004  meshPointsOnThisSegment,
1005  facePoint
1006  )
1007  ==
1008  -1
1009  )
1010  {
1011  meshPointsOnThisSegment.append(facePoint);
1012  }
1013  }
1014  }
1015 
1016  meshEdgesOnThisSegment.shrink();
1017 
1018  meshPointsOnThisSegment.shrink();
1019 
1020  if (meshFacesOnThisSegment.size())
1021  {
1022  if (faceT == -1)
1023  {
1024  FatalErrorIn ("referredCellList.C")
1025  << nl << "faceT == -1 encountered but "
1026  << meshFacesOnThisSegment.size()
1027  << " faces found on patch segment."
1028  << abort(FatalError);
1029  }
1030 
1031  if (iterationNo == 0)
1032  {
1033  // Assessing real cells in range is only required on
1034  // the 1st iteration because they do not change from
1035  // iteration to iteration.
1036 
1037  labelList realCellsFoundInRange
1038  (
1040  (
1041  meshFacesOnThisSegment,
1042  meshEdgesOnThisSegment,
1043  meshPointsOnThisSegment
1044  )
1045  );
1046 
1047  forAll(realCellsFoundInRange,cFIR)
1048  {
1049  const label realCell =
1050  realCellsFoundInRange[cFIR];
1051 
1052  referredCell cellToRefer
1053  (
1054  mesh,
1056  realCell,
1057  patch.faceCentres()[faceT],
1058  neighbFaceCentres[faceT],
1059  patch.faceNormals()[faceT],
1060  neighbFaceAreas[faceT]
1061  /(mag(neighbFaceAreas[faceT]) + VSMALL)
1062  );
1063 
1064  referredCellsToTransfer.append(cellToRefer);
1065 
1066  // add real cells found in range of processor
1067  // patch to whole mesh list
1068 
1069  if (findIndex (rCellsWRRP, realCell) == -1)
1070  {
1071  rCellsWRRP.append(realCell);
1072  }
1073  }
1074  }
1075 
1076  referredInteractionList.shrink();
1077 
1078  labelList referredCellsFoundInRange
1079  (
1081  (
1082  referredInteractionList,
1083  meshFacesOnThisSegment,
1084  meshEdgesOnThisSegment,
1085  meshPointsOnThisSegment
1086  )
1087  );
1088 
1089  forAll(referredCellsFoundInRange,cFIR)
1090  {
1091  referredCell& existingRefCell =
1092  referredInteractionList
1093  [
1094  referredCellsFoundInRange[cFIR]
1095  ];
1096 
1097  referredCell cellToReRefer =
1098  existingRefCell.reRefer
1099  (
1100  patch.faceCentres()[faceT],
1101  neighbFaceCentres[faceT],
1102  patch.faceNormals()[faceT],
1103  neighbFaceAreas[faceT]
1104  /(mag(neighbFaceAreas[faceT]) + VSMALL)
1105  );
1106 
1107  referredCellsToTransfer.append(cellToReRefer);
1108  }
1109  }
1110  }
1111 
1112  referredCellsToTransfer.shrink();
1113 
1114  // Send these cells to the neighbouring processor.
1115 
1116  {
1117  OPstream toNeighbProc
1118  (
1120  patch.neighbProcNo()
1121  );
1122 
1123  toNeighbProc << referredCellsToTransfer;
1124  }
1125  }
1126 
1127  forAll(procPatches,pP)
1128  {
1129  const processorPolyPatch& patch =
1130  refCast<const processorPolyPatch>
1131  (
1132  mesh.boundaryMesh()[procPatches[pP]]
1133  );
1134 
1135  // Receive the cells from neighbour
1136 
1137  List<referredCell> referredCellsFromNeighbour(patch.size());
1138 
1139  {
1140  IPstream fromNeighbProc
1141  (
1143  patch.neighbProcNo()
1144  );
1145 
1146  fromNeighbProc >> referredCellsFromNeighbour;
1147  }
1148 
1149  // Check to see if they are duplicates, if not append
1150  // them to the referredInteractionList
1151 
1152  forAll(referredCellsFromNeighbour,rCFN)
1153  {
1154  referredCell& cellToRefer =
1155  referredCellsFromNeighbour[rCFN];
1156 
1157  // Test all existing referred and real cells to check
1158  // duplicates are not being made or cells aren't being
1159  // referred back onto themselves
1160 
1161  bool addCellToRefer = true;
1162 
1163  // Check if cellToRefer is an existing referred cell
1164 
1165  forAll(referredInteractionList, rIL)
1166  {
1167  if (cellToRefer.duplicate(referredInteractionList[rIL]))
1168  {
1169  addCellToRefer = false;
1170 
1171  break;
1172  }
1173  }
1174 
1175  // Check for cellToRefer being referred back ontop of a real
1176  // cell
1177 
1178  if
1179  (
1180  cellToRefer.duplicate
1181  (
1183  mesh.nCells()
1184  )
1185  )
1186  {
1187  addCellToRefer = false;
1188  }
1189 
1190  if (addCellToRefer)
1191  {
1192  referredInteractionList.append(cellToRefer);
1193  }
1194  }
1195  }
1196  }
1197 
1198  if (iterationNo == 0)
1199  {
1200  // record all real cells in range of any referring patch (cyclic or
1201  // processor) on the first iteration when the real cells are
1202  // evaluated.
1203 
1204  rCellsWRRP.shrink();
1205 
1206  // construct {faces, edges, points}WithinRCutMaxOfAnyReferringPatch
1207 
1208  forAll(rCellsWRRP, rCWR)
1209  {
1210  const label realCell(rCellsWRRP[rCWR]);
1211 
1212  const labelList& rCFaces
1213  (
1214  mesh.cells()[realCell]
1215  );
1216 
1217  forAll(rCFaces, rCF)
1218  {
1219  const label f(rCFaces[rCF]);
1220 
1221  if (findIndex(rFacesWRRP,f) == -1)
1222  {
1223  rFacesWRRP.append(f);
1224  }
1225  }
1226 
1227  const labelList& rCEdges
1228  (
1229  mesh.cellEdges()[realCell]
1230  );
1231 
1232  forAll(rCEdges, rCE)
1233  {
1234  const label e(rCEdges[rCE]);
1235 
1236  if (findIndex(rEdgesWRRP,e) == -1)
1237  {
1238  rEdgesWRRP.append(e);
1239  }
1240  }
1241 
1242  const labelList& rCPoints
1243  (
1244  mesh.cellPoints()[realCell]
1245  );
1246 
1247  forAll(rCPoints, rCP)
1248  {
1249  const label p(rCPoints[rCP]);
1250 
1251  if (findIndex(rPointsWRRP,p) == -1)
1252  {
1253  rPointsWRRP.append(p);
1254  }
1255  }
1256  }
1257 
1258  rFacesWRRP.shrink();
1259 
1260  rEdgesWRRP.shrink();
1261 
1262  rPointsWRRP.shrink();
1263  }
1264 
1265  iterationNo++;
1266 
1267  cellsReferredThisIteration =
1268  referredInteractionList.size() - refIntListStartSize;
1269 
1270  reduce(cellsReferredThisIteration, sumOp<label>());
1271 
1272  Info<< tab << "Cells added this iteration: "
1273  << cellsReferredThisIteration << endl;
1274  }
1275 
1276  referredInteractionList.shrink();
1277 
1278  // Info<< "referredInteractionList.size() = "
1279  // << referredInteractionList.size() << endl;
1280 
1281  // forAll(referredInteractionList,rIL)
1282  // {
1283  // Info<< referredInteractionList[rIL];
1284  // }
1285 
1286  (*this).setSize
1287  (
1288  referredInteractionList.size()
1289  );
1290 
1291  forAll(referredInteractionList, rIL)
1292  {
1293  (*this)[rIL] = referredInteractionList[rIL];
1294  }
1295 
1296  Info<< nl << "Finding real cells in range of referred cells" << endl;
1297 
1298  forAll(*this, rC)
1299  {
1300  const polyMesh& mesh(il_.mesh());
1301 
1302  referredCell& refCell = (*this)[rC];
1303 
1304  DynamicList<label> realCellsFoundInRange;
1305 
1306  const vectorList& refCellPoints = refCell.vertexPositions();
1307 
1308  forAll(rFacesWRRP, rCF)
1309  {
1310  const label f(rFacesWRRP[rCF]);
1311 
1312  if (il_.testPointFaceDistance(refCellPoints,f))
1313  {
1314  const label cellO(mesh.faceOwner()[f]);
1315 
1316  if (findIndex(realCellsFoundInRange, cellO) == -1)
1317  {
1318  realCellsFoundInRange.append(cellO);
1319  }
1320 
1321  if (mesh.isInternalFace(f))
1322  {
1323  // boundary faces will not have neighbour information
1324 
1325  const label cellN(mesh.faceNeighbour()[f]);
1326 
1327  if (findIndex(realCellsFoundInRange, cellN) == -1)
1328  {
1329  realCellsFoundInRange.append(cellN);
1330  }
1331  }
1332  }
1333  }
1334 
1335  forAll(rPointsWRRP, rCP)
1336  {
1337  const label p(rPointsWRRP[rCP]);
1338 
1339  if (il_.testPointFaceDistance(p,refCell))
1340  {
1341  const labelList& pCells(mesh.pointCells()[p]);
1342 
1343  forAll(pCells, pC)
1344  {
1345  const label cellI(pCells[pC]);
1346 
1347  if (findIndex(realCellsFoundInRange, cellI) == -1)
1348  {
1349  realCellsFoundInRange.append(cellI);
1350  }
1351  }
1352  }
1353  }
1354 
1355 
1356  const edgeList& refCellEdges = refCell.edges();
1357 
1358  forAll(rEdgesWRRP, rCE)
1359  {
1360  const label edgeIIndex(rEdgesWRRP[rCE]);
1361 
1362  const edge& eI(mesh.edges()[edgeIIndex]);
1363 
1364  forAll(refCellEdges, rCE)
1365  {
1366  const edge& eJ(refCellEdges[rCE]);
1367 
1368  if
1369  (
1371  (
1372  eI,
1373  refCellPoints[eJ.start()],
1374  refCellPoints[eJ.end()]
1375  )
1376  )
1377  {
1378  const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
1379 
1380  forAll(eICells, eIC)
1381  {
1382  const label cellI(eICells[eIC]);
1383 
1384  if (findIndex(realCellsFoundInRange, cellI) == -1)
1385  {
1386  realCellsFoundInRange.append(cellI);
1387  }
1388  }
1389  }
1390  }
1391  }
1392 
1393 // scalar rCutMaxSqr = molCloud_.rCutMax()*molCloud_.rCutMax();
1394 //
1395 // forAll (molCloud_.mesh().points(), pointIIndex)
1396 // {
1397 // const point& ptI
1398 // (
1399 // molCloud_.mesh().points()[pointIIndex]
1400 // );
1401 //
1402 // forAll(refCellPoints, rCP)
1403 // {
1404 // if (magSqr(ptI - refCellPoints[rCP]) <= rCutMaxSqr)
1405 // {
1406 // const labelList& ptICells
1407 // (
1408 // molCloud_.mesh().pointCells()[pointIIndex]
1409 // );
1410 //
1411 // forAll(ptICells, pIC)
1412 // {
1413 // const label cellI(ptICells[pIC]);
1414 //
1415 // if (findIndex(realCellsFoundInRange, cellI) == -1)
1416 // {
1417 // realCellsFoundInRange.append(cellI);
1418 // }
1419 // }
1420 // }
1421 // }
1422 // }
1423 
1424  refCell.realCells() = realCellsFoundInRange.shrink();
1425  }
1426 }
1427 
1428 
1429 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1430 
1433  interactionLists& il,
1434  bool pointPointListBuild
1435 )
1436 :
1438  il_(il)
1439 {
1440  buildReferredCellList(pointPointListBuild);
1441 }
1442 
1443 
1445 :
1446  List<referredCell>(),
1447  il_(il)
1448 {
1449  Info<< "Read referredCellList from disk not implemented" << endl;
1450 }
1451 
1452 
1453 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1454 
1456 {}
1457 
1458 
1459 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1460 
1464 )
1465 {
1466  // Create referred molecules for sending using cell occupancy and
1467  // cellSendingReferralLists
1468 
1469  forAll(il_.cellSendingReferralLists(), cSRL)
1470  {
1471  const sendingReferralList& sRL
1472  (
1473  il_.cellSendingReferralLists()[cSRL]
1474  );
1475 
1476  List<DynamicList<referredMolecule> > molsToReferOut(sRL.size());
1477 
1478  forAll(sRL, sRLI)
1479  {
1480  List<molecule*> realMols = cellOccupancy[sRL[sRLI]];
1481 
1482  forAll (realMols, rM)
1483  {
1484  molecule* mol = realMols[rM];
1485 
1486  molsToReferOut[sRLI].append
1487  (
1489  (
1490  mol->id(),
1491  mol->position(),
1492  mol->sitePositions()
1493  )
1494  );
1495  }
1496 
1497  molsToReferOut[sRLI].shrink();
1498  }
1499 
1500  // Send lists of referred molecules to other processors
1501 
1502  if (sRL.destinationProc() != Pstream::myProcNo())
1503  {
1504  if (Pstream::parRun())
1505  {
1506  OPstream toInteractingProc
1507  (
1509  sRL.destinationProc()
1510  );
1511 
1512  toInteractingProc << molsToReferOut;
1513  }
1514  }
1515  else
1516  {
1517  // Refer molecules directly for referred cells on the same
1518  // processor.
1519 
1520  const receivingReferralList& rRL
1521  (
1522  il_.cellReceivingReferralLists()[cSRL]
1523  );
1524 
1525  forAll(rRL, rRLI)
1526  {
1527  forAll(rRL[rRLI], rC)
1528  {
1529  referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
1530 
1531  refCellToRefMolsTo.referInMols(molsToReferOut[rRLI]);
1532  }
1533  }
1534  }
1535  }
1536 
1537  // Receive referred molecule lists to and distribute to referredCells
1538  // according tocellReceivingReferralLists, referredCells deal with the
1539  // transformations themselves
1540 
1541  forAll(il_.cellReceivingReferralLists(), cRRL)
1542  {
1543  const receivingReferralList& rRL
1544  (
1545  il_.cellReceivingReferralLists()[cRRL]
1546  );
1547 
1548  List<List<referredMolecule> > molsToReferIn(rRL.size());
1549 
1550  if (rRL.sourceProc() != Pstream::myProcNo())
1551  {
1552  if (Pstream::parRun())
1553  {
1554  IPstream fromInteractingProc
1555  (
1557  rRL.sourceProc()
1558  );
1559 
1560  fromInteractingProc >> molsToReferIn;
1561  }
1562 
1563  forAll(rRL, rRLI)
1564  {
1565  forAll(rRL[rRLI], rC)
1566  {
1567  referredCell& refCellToRefMolsTo = (*this)[rRL[rRLI][rC]];
1568 
1569  refCellToRefMolsTo.referInMols(molsToReferIn[rRLI]);
1570  }
1571  }
1572  }
1573  }
1574 }
1575 
1576 
1577 // ************************ vim: set sw=4 sts=4 et: ************************ //