FreeFOAM The Cross-Platform CFD Toolkit
interactionLists.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 "interactionLists.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 Foam::scalar Foam::interactionLists::transTol = 1e-12;
31 
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 void Foam::interactionLists::buildCellReferralLists()
36 {
37  Info<< nl << "Determining molecule referring schedule" << endl;
38 
39  const referredCellList& refIntL(ril());
40 
41  DynamicList<label> referralProcs;
42 
43  // Run through all referredCells to build list of interacting processors
44 
45  forAll(refIntL, rIL)
46  {
47  const referredCell& rC(refIntL[rIL]);
48 
49  if (findIndex(referralProcs, rC.sourceProc()) == -1)
50  {
51  referralProcs.append(rC.sourceProc());
52  }
53  }
54 
55  referralProcs.shrink();
56 
57 // Pout << "referralProcs: " << nl << referralProcs << endl;
58 
59  List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
60 
61  // Use autoPtr to help GCC compilers (< 4.3) which suffer from core
62  // language defect 391, i.e. require a copy-ctor to be present when passing
63  // a temporary as an rvalue.
64  List<DynamicList<autoPtr<DynamicList<label> > > >
65  cellReceivingReferralLists(referralProcs.size());
66 
67  // Run through all referredCells again building up send and receive info
68 
69  forAll(refIntL, rIL)
70  {
71  const referredCell& rC(refIntL[rIL]);
72 
73  label rPI = findIndex(referralProcs, rC.sourceProc());
74 
75  DynamicList<autoPtr<DynamicList<label> > >& rRL(cellReceivingReferralLists[rPI]);
76 
77  DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
78 
79  label existingSource = findIndex(sRL, rC.sourceCell());
80 
81  // Check to see if this source cell has already been allocated to
82  // come to this processor. If not, add the source cell to the sending
83  // list and add the current referred cell to the receiving list.
84 
85  // It shouldn't be possible for the sending and receiving lists to be
86  // different lengths, because their append operations happen at the
87  // same time.
88 
89  if (existingSource == -1)
90  {
91  sRL.append(rC.sourceCell());
92 
93  rRL.append
94  (
95  autoPtr<DynamicList<label> >(new DynamicList<label> (labelList(1,rIL)))
96  );
97  }
98  else
99  {
100  rRL[existingSource]->append(rIL);
101 
102  rRL[existingSource]->shrink();
103  }
104  }
105 
106  forAll(referralProcs, rPI)
107  {
108  DynamicList<autoPtr<DynamicList<label> > >& rRL(cellReceivingReferralLists[rPI]);
109 
110  DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
111 
112  sRL.shrink();
113 
114  rRL.shrink();
115  }
116 
117  // It is assumed that cell exchange is reciprocal, if proc A has cells to
118  // send to proc B, then proc B must have some to send to proc A.
119 
120  cellReceivingReferralLists_.setSize(referralProcs.size());
121 
122  cellSendingReferralLists_.setSize(referralProcs.size());
123 
124  forAll(referralProcs, rPI)
125  {
126  DynamicList<autoPtr<DynamicList<label> > >& rRL(cellReceivingReferralLists[rPI]);
127 
128  labelListList translLL(rRL.size());
129 
130  forAll(rRL, rRLI)
131  {
132  translLL[rRLI] = rRL[rRLI]();
133  }
134 
135  cellReceivingReferralLists_[rPI] = receivingReferralList
136  (
137  referralProcs[rPI],
138  translLL
139  );
140  }
141 
142  // Send sendingReferralLists to each interacting processor.
143 
144  forAll(referralProcs, rPI)
145  {
146 
147  DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
148 
149  if (referralProcs[rPI] != Pstream::myProcNo())
150  {
151  if (Pstream::parRun())
152  {
153  OPstream toInteractingProc
154  (
156  referralProcs[rPI]
157  );
158 
159  toInteractingProc << sendingReferralList
160  (
162  sRL
163  );
164  }
165  }
166  }
167 
168  // Receive sendingReferralLists from each interacting processor.
169 
170  forAll(referralProcs, rPI)
171  {
172  if (referralProcs[rPI] != Pstream::myProcNo())
173  {
174  if (Pstream::parRun())
175  {
176  IPstream fromInteractingProc
177  (
179  referralProcs[rPI]
180  );
181 
182  fromInteractingProc >> cellSendingReferralLists_[rPI];
183  }
184  }
185  else
186  {
187  cellSendingReferralLists_[rPI] = sendingReferralList
188  (
191  );
192  }
193  }
194 }
195 
196 
197 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
198 
199 
200 Foam::interactionLists::interactionLists
201 (
202  const polyMesh& mesh,
203  scalar rCutMaxSqr,
204  bool pointPointListBuild
205 )
206 :
207  mesh_(mesh),
208  rCutMaxSqr_(rCutMaxSqr),
209  dil_(*this, pointPointListBuild),
210  ril_(*this, pointPointListBuild),
211  cellSendingReferralLists_(),
212  cellReceivingReferralLists_()
213 {
214  buildCellReferralLists();
215 }
216 
217 
218 Foam::interactionLists::interactionLists(const polyMesh& mesh)
219 :
220  mesh_(mesh),
221  dil_(*this),
222  ril_(*this)
223 {}
224 
225 
226 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
227 
229 {}
230 
231 
232 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
233 
235 (
236  const label ptI,
237  const label ptJ
238 ) const
239 {
240  return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_);
241 }
242 
243 
245 (
246  const edge& eI,
247  const edge& eJ
248 ) const
249 {
250  const vector& eJs(mesh_.points()[eJ.start()]);
251  const vector& eJe(mesh_.points()[eJ.end()]);
252 
253  return testEdgeEdgeDistance(eI, eJs, eJe);
254 }
255 
256 
258 (
259  const label p,
260  const label faceNo
261 ) const
262 {
263  const vector& pointPosition(mesh_.points()[p]);
264 
265  return testPointFaceDistance(pointPosition, faceNo);
266 }
267 
268 
270 (
271  const label p,
272  const referredCell& refCell
273 ) const
274 {
275  const vector& pointPosition(mesh_.points()[p]);
276 
277  forAll (refCell.faces(), rCF)
278  {
279  if
280  (
281  testPointFaceDistance
282  (
283  pointPosition,
284  refCell.faces()[rCF],
285  refCell.vertexPositions(),
286  refCell.faceCentres()[rCF],
287  refCell.faceAreas()[rCF]
288  )
289  )
290  {
291  return true;
292  }
293  }
294 
295  return false;
296 }
297 
298 
300 (
301  const vectorList& pointsToTest,
302  const label faceNo
303 ) const
304 {
305  forAll(pointsToTest, pTT)
306  {
307  const vector& p(pointsToTest[pTT]);
308 
309  // if any point in the list is in range of the face
310  // then the rest do not need to be tested and
311  // true can be returned
312 
313  if (testPointFaceDistance(p, faceNo))
314  {
315  return true;
316  }
317  }
318 
319  return false;
320 }
321 
322 
324 (
325  const vector& p,
326  const label faceNo
327 ) const
328 {
329  const face& faceToTest(mesh_.faces()[faceNo]);
330 
331  const vector& faceC(mesh_.faceCentres()[faceNo]);
332 
333  const vector& faceA(mesh_.faceAreas()[faceNo]);
334 
335  const vectorList& points(mesh_.points());
336 
337  return testPointFaceDistance
338  (
339  p,
340  faceToTest,
341  points,
342  faceC,
343  faceA
344  );
345 }
346 
347 
349 (
350  const vector& p,
351  const labelList& faceToTest,
352  const vectorList& points,
353  const vector& faceC,
354  const vector& faceA
355 ) const
356 {
357  vector faceN(faceA/mag(faceA));
358 
359  scalar perpDist((p - faceC) & faceN);
360 
361  if (magSqr(perpDist) > rCutMaxSqr_)
362  {
363  return false;
364  }
365 
366  vector pointOnPlane = (p - faceN * perpDist);
367 
368  if (magSqr(faceC - pointOnPlane) < rCutMaxSqr_*1e-8)
369  {
370  // If pointOnPlane is very close to the face centre
371  // then defining the local axes will be inaccurate
372  // and it is very likely that pointOnPlane will be
373  // inside the face, so return true if the points
374  // are in range to be safe
375 
376  return (magSqr(pointOnPlane - p) <= rCutMaxSqr_);
377  }
378 
379  vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
380 
381  vector yAxis =
382  ((faceC - pointOnPlane) ^ faceN)
383  /mag((faceC - pointOnPlane) ^ faceN);
384 
385  List<vector2D> local2DVertices(faceToTest.size());
386 
387  forAll(faceToTest, fTT)
388  {
389  const vector& V(points[faceToTest[fTT]]);
390 
391  if (magSqr(V-p) <= rCutMaxSqr_)
392  {
393  return true;
394  }
395 
396  local2DVertices[fTT] = vector2D
397  (
398  ((V - pointOnPlane) & xAxis),
399  ((V - pointOnPlane) & yAxis)
400  );
401  }
402 
403  scalar localFaceCx((faceC - pointOnPlane) & xAxis);
404 
405  scalar la_valid = -1;
406 
407  forAll(local2DVertices, fV)
408  {
409  const vector2D& va(local2DVertices[fV]);
410 
411  const vector2D& vb
412  (
413  local2DVertices[(fV + 1) % local2DVertices.size()]
414  );
415 
416  if (mag(vb.y()-va.y()) > SMALL)
417  {
418  scalar la =
419  (
420  va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
421  )
422  /localFaceCx;
423 
424  scalar lv = -va.y()/(vb.y() - va.y());
425 
426 
427  if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
428  {
429  la_valid = la;
430 
431  break;
432  }
433  }
434  }
435 
436  if (la_valid < 0)
437  {
438  // perpendicular point inside face, nearest point is pointOnPlane;
439  return (magSqr(pointOnPlane-p) <= rCutMaxSqr_);
440  }
441  else
442  {
443  // perpendicular point outside face, nearest point is
444  // on edge that generated la_valid;
445  return
446  (
447  magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
448  <= rCutMaxSqr_
449  );
450  }
451 
452  // if the algorithm hasn't returned anything by now then something has
453  // gone wrong.
454 
455  FatalErrorIn("interactionLists.C") << nl
456  << "point " << p << " to face " << faceToTest
457  << " comparison did not find a nearest point"
458  << " to be inside or outside face."
459  << abort(FatalError);
460 
461  return false;
462 }
463 
464 
466 (
467  const edge& eI,
468  const vector& eJs,
469  const vector& eJe
470 ) const
471 {
472  vector a(eI.vec(mesh_.points()));
473  vector b(eJe - eJs);
474 
475  const vector& eIs(mesh_.points()[eI.start()]);
476 
477  vector c(eJs - eIs);
478 
479  vector crossab = a ^ b;
480  scalar magCrossSqr = magSqr(crossab);
481 
482  if (magCrossSqr > VSMALL)
483  {
484  // If the edges are parallel then a point-face
485  // search will pick them up
486 
487  scalar s = ((c ^ b) & crossab)/magCrossSqr;
488  scalar t = ((c ^ a) & crossab)/magCrossSqr;
489 
490  // Check for end points outside of range 0..1
491  // If the closest point is outside this range
492  // a point-face search will have found it.
493 
494  return
495  (
496  s >= 0
497  && s <= 1
498  && t >= 0
499  && t <= 1
500  && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr_
501  );
502  }
503 
504  return false;
505 }
506 
507 
509 (
510  const labelList& segmentFaces,
511  const labelList& segmentEdges,
512  const labelList& segmentPoints
513 ) const
514 {
515  DynamicList<label> realCellsFoundInRange;
516 
517  forAll(segmentFaces, sF)
518  {
519  const label f = segmentFaces[sF];
520 
521  forAll (mesh_.points(), p)
522  {
523  if (testPointFaceDistance(p, f))
524  {
525  const labelList& pCells(mesh_.pointCells()[p]);
526 
527  forAll(pCells, pC)
528  {
529  const label cellI(pCells[pC]);
530 
531  if (findIndex(realCellsFoundInRange, cellI) == -1)
532  {
533  realCellsFoundInRange.append(cellI);
534  }
535  }
536  }
537  }
538  }
539 
540  forAll(segmentPoints, sP)
541  {
542  const label p = segmentPoints[sP];
543 
544  forAll(mesh_.faces(), f)
545  {
546  if (testPointFaceDistance(p, f))
547  {
548  const label cellO(mesh_.faceOwner()[f]);
549 
550  if (findIndex(realCellsFoundInRange, cellO) == -1)
551  {
552  realCellsFoundInRange.append(cellO);
553  }
554 
555  if (mesh_.isInternalFace(f))
556  {
557  // boundary faces will not have neighbour information
558 
559  const label cellN(mesh_.faceNeighbour()[f]);
560 
561  if (findIndex(realCellsFoundInRange, cellN) == -1)
562  {
563  realCellsFoundInRange.append(cellN);
564  }
565  }
566  }
567  }
568  }
569 
570  forAll(segmentEdges, sE)
571  {
572  const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
573 
574  forAll (mesh_.edges(), edgeIIndex)
575  {
576  const edge& eI(mesh_.edges()[edgeIIndex]);
577 
578  if (testEdgeEdgeDistance(eI, eJ))
579  {
580  const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
581 
582  forAll(eICells, eIC)
583  {
584  const label cellI(eICells[eIC]);
585 
586  if (findIndex(realCellsFoundInRange, cellI) == -1)
587  {
588  realCellsFoundInRange.append(cellI);
589  }
590  }
591  }
592  }
593  }
594 
595  return realCellsFoundInRange.shrink();
596 }
597 
598 
600 (
601  const List<referredCell>& referredInteractionList,
602  const labelList& segmentFaces,
603  const labelList& segmentEdges,
604  const labelList& segmentPoints
605 ) const
606 {
607  DynamicList<label> referredCellsFoundInRange;
608 
609  forAll(segmentFaces, sF)
610  {
611  const label f = segmentFaces[sF];
612 
613  forAll(referredInteractionList, rIL)
614  {
615  const vectorList& refCellPoints
616  = referredInteractionList[rIL].vertexPositions();
617 
618  if (testPointFaceDistance(refCellPoints, f))
619  {
620  if (findIndex(referredCellsFoundInRange, rIL) == -1)
621  {
622  referredCellsFoundInRange.append(rIL);
623  }
624  }
625  }
626  }
627 
628  forAll(segmentPoints, sP)
629  {
630  const label p = segmentPoints[sP];
631 
632  forAll(referredInteractionList, rIL)
633  {
634  const referredCell& refCell(referredInteractionList[rIL]);
635 
636  if (testPointFaceDistance(p, refCell))
637  {
638  if (findIndex(referredCellsFoundInRange, rIL) == -1)
639  {
640  referredCellsFoundInRange.append(rIL);
641  }
642  }
643  }
644  }
645 
646  forAll(segmentEdges, sE)
647  {
648  const edge& eI(mesh_.edges()[segmentEdges[sE]]);
649 
650  forAll(referredInteractionList, rIL)
651  {
652  const vectorList& refCellPoints
653  = referredInteractionList[rIL].vertexPositions();
654 
655  const edgeList& refCellEdges
656  = referredInteractionList[rIL].edges();
657 
658  forAll(refCellEdges, rCE)
659  {
660  const edge& eJ(refCellEdges[rCE]);
661 
662  if
663  (
664  testEdgeEdgeDistance
665  (
666  eI,
667  refCellPoints[eJ.start()],
668  refCellPoints[eJ.end()]
669  )
670  )
671  {
672  if(findIndex(referredCellsFoundInRange, rIL) == -1)
673  {
674  referredCellsFoundInRange.append(rIL);
675  }
676  }
677  }
678  }
679  }
680 
681  return referredCellsFoundInRange.shrink();
682 }
683 
684 
685 // ************************ vim: set sw=4 sts=4 et: ************************ //