rayShooting.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) 2013-2015 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 "rayShooting.H"
28 #include "triSurfaceMesh.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 defineTypeNameAndDebug(rayShooting, 0);
38 addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 void rayShooting::splitLine
44 (
45  const line<point, point>& l,
46  const scalar& pert,
47  DynamicList<Vb::Point>& initialPoints
48 ) const
49 {
50  Foam::point midPoint(l.centre());
51  const scalar localCellSize(cellShapeControls().cellSize(midPoint));
52 
53  const scalar minDistFromSurfaceSqr
54  (
56  *sqr(localCellSize)
57  );
58 
59  if
60  (
61  magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
62  && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
63  )
64  {
65  // Add extra points if line length is much bigger than local cell size
66 // const scalar lineLength(l.mag());
67 //
68 // if (lineLength > 4.0*localCellSize)
69 // {
70 // splitLine
71 // (
72 // line<point, point>(l.start(), midPoint),
73 // pert,
74 // initialPoints
75 // );
76 //
77 // splitLine
78 // (
79 // line<point, point>(midPoint, l.end()),
80 // pert,
81 // initialPoints
82 // );
83 // }
84 
85  if (randomiseInitialGrid_)
86  {
87  Foam::point newPt
88  (
89  midPoint.x() + pert*(rndGen().scalar01() - 0.5),
90  midPoint.y() + pert*(rndGen().scalar01() - 0.5),
91  midPoint.z() + pert*(rndGen().scalar01() - 0.5)
92  );
93 
94  if
95  (
97  (
98  midPoint,
99  newPt
100  )
101  )
102  {
103  midPoint = newPt;
104  }
105  else
106  {
108  << "Point perturbation crosses a surface. Not inserting."
109  << endl;
110  }
111  }
112 
113  initialPoints.append(toPoint(midPoint));
114  }
115 }
116 
117 
118 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
119 
121 (
122  const dictionary& initialPointsDict,
123  const Time& runTime,
124  Random& rndGen,
125  const conformationSurfaces& geometryToConformTo,
126  const cellShapeControl& cellShapeControls,
127  const autoPtr<backgroundMeshDecomposition>& decomposition
128 )
129 :
131  (
132  typeName,
133  initialPointsDict,
134  runTime,
135  rndGen,
136  geometryToConformTo,
137  cellShapeControls,
138  decomposition
139  ),
140  randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
141  randomPerturbationCoeff_
142  (
143  readScalar(detailsDict().lookup("randomPerturbationCoeff"))
144  )
145 {}
146 
147 
148 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
149 
150 List<Vb::Point> rayShooting::initialPoints() const
151 {
152  // Loop over surface faces
153  const searchableSurfaces& surfaces = geometryToConformTo().geometry();
154  const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
155 
156  const scalar maxRayLength(surfaces.bounds().mag());
157 
158  // Initialise points list
159  label initialPointsSize = 0;
160  forAll(surfaces, surfI)
161  {
162  initialPointsSize += surfaces[surfI].size();
163  }
164 
165  DynamicList<Vb::Point> initialPoints(initialPointsSize);
166 
167  forAll(surfacesToConformTo, surfI)
168  {
169  const searchableSurface& s = surfaces[surfacesToConformTo[surfI]];
170 
171  tmp<pointField> faceCentresTmp(s.coordinates());
172  const pointField& faceCentres = faceCentresTmp();
173 
174  Info<< " Shoot rays from " << s.name() << nl
175  << " nRays = " << faceCentres.size() << endl;
176 
177  forAll(faceCentres, fcI)
178  {
179  const Foam::point& fC = faceCentres[fcI];
180 
181  if
182  (
184  && !decomposition().positionOnThisProcessor(fC)
185  )
186  {
187  continue;
188  }
189 
190  const scalar pert
191  (
192  randomPerturbationCoeff_
193  *cellShapeControls().cellSize(fC)
194  );
195 
196  pointIndexHit surfHitStart;
197  label hitSurfaceStart;
198 
199  // Face centres should be on the surface so search distance can be
200  // small
202  (
203  fC,
204  sqr(pert),
205  surfHitStart,
206  hitSurfaceStart
207  );
208 
209  vectorField normStart(1, vector::min);
211  (
212  hitSurfaceStart,
213  List<pointIndexHit>(1, surfHitStart),
214  normStart
215  );
216 
217  pointIndexHit surfHitEnd;
218  label hitSurfaceEnd;
219 
221  (
222  fC - normStart[0]*pert,
223  fC - normStart[0]*maxRayLength,
224  surfHitEnd,
225  hitSurfaceEnd
226  );
227 
228  if (surfHitEnd.hit())
229  {
230  vectorField normEnd(1, vector::min);
232  (
233  hitSurfaceEnd,
234  List<pointIndexHit>(1, surfHitEnd),
235  normEnd
236  );
237 
238  if ((normStart[0] & normEnd[0]) < 0)
239  {
240  line<point, point> l(fC, surfHitEnd.hitPoint());
241 
242  if (Pstream::parRun())
243  {
244  // Clip the line in parallel
245  pointIndexHit procIntersection =
247  (
248  l.start(),
249  l.end()
250  );
251 
252  if (procIntersection.hit())
253  {
254  l =
255  line<point, point>
256  (
257  l.start(),
258  procIntersection.hitPoint()
259  );
260  }
261  }
262 
263  splitLine
264  (
265  l,
266  pert,
267  initialPoints
268  );
269  }
270  }
271  }
272  }
273 
274  return initialPoints.shrink();
275 }
276 
277 
278 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
279 
280 } // End namespace Foam
281 
282 // ************************************************************************* //
virtual List< Vb::Point > initialPoints() const
Return the initial points for the conformalVoronoiMesh.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const cellShapeControl & cellShapeControls() const
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
virtual List< Vb::Point > initialPoints() const
Return the initial points for the conformalVoronoiMesh.
rayShooting(const dictionary &initialPointsDict, const Time &runTime, Random &rndGen, const conformationSurfaces &geometryToConformTo, const cellShapeControl &cellShapeControls, const autoPtr< backgroundMeshDecomposition > &decomposition)
Construct from components.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
void getNormal(const label hitSurface, const List< pointIndexHit > &surfHit, vectorField &normal) const
const dictionary & detailsDict() const
Const access to the details dictionary.
Macros for easy insertion into run-time selection tables.
dictionary()
Construct top-level dictionary null.
Definition: dictionary.C:113
bool findSurfaceAnyIntersection(const point &start, const point &end) const
const labelList & surfaces() const
Return the surface indices.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
initialPointsMethod(const word &type, const dictionary &initialPointsDict, const Time &runTime, Random &rndGen, const conformationSurfaces &geometryToConformTo, const cellShapeControl &cellShapeControls, const autoPtr< backgroundMeshDecomposition > &decomposition)
Construct from components.
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
PointFrompoint toPoint(const Foam::point &p)
List< label > labelList
A List of labels.
Definition: labelList.H:56
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
void findSurfaceNearestIntersection(const point &start, const point &end, pointIndexHit &surfHit, label &hitSurface) const
Finding the nearestIntersection of the surface to start.
#define WarningInFunction
Report a warning using Foam::Warning.
const conformationSurfaces & geometryToConformTo() const
messageStream Info
Field< vector > vectorField
Specialisation of Field<T> for vector.
const backgroundMeshDecomposition & decomposition() const
scalar scalar01()
Scalar [0..1] (so including 0,1)
Definition: Random.C:67
scalar minimumSurfaceDistanceCoeffSqr_
Only allow the placement of initial points that are within the.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451