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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2015 OpenFOAM Foundation
9 Copyright (C) 2018 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "rayShooting.H"
31#include "triSurfaceMesh.H"
32
33// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34
35namespace Foam
36{
37 defineTypeNameAndDebug(rayShooting, 0);
38 addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
39}
40
41
42// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43
44void Foam::rayShooting::splitLine
45(
46 const line<point, point>& l,
47 const scalar& pert,
48 DynamicList<Vb::Point>& initialPoints
49) const
50{
51 Foam::point midPoint(l.centre());
52 const scalar localCellSize(cellShapeControls().cellSize(midPoint));
53
54 const scalar minDistFromSurfaceSqr
55 (
57 *sqr(localCellSize)
58 );
59
60 if
61 (
62 magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
63 && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
64 )
65 {
66 // Add extra points if line length is much bigger than local cell size
67// const scalar lineLength(l.mag());
68//
69// if (lineLength > 4.0*localCellSize)
70// {
71// splitLine
72// (
73// line<point, point>(l.start(), midPoint),
74// pert,
75// initialPoints
76// );
77//
78// splitLine
79// (
80// line<point, point>(midPoint, l.end()),
81// pert,
82// initialPoints
83// );
84// }
85
86 if (randomiseInitialGrid_)
87 {
88 Foam::point newPt
89 (
90 midPoint.x() + pert*(rndGen().sample01<scalar>() - 0.5),
91 midPoint.y() + pert*(rndGen().sample01<scalar>() - 0.5),
92 midPoint.z() + pert*(rndGen().sample01<scalar>() - 0.5)
93 );
94
95 if
96 (
97 !geometryToConformTo().findSurfaceAnyIntersection
98 (
99 midPoint,
100 newPt
101 )
102 )
103 {
104 midPoint = newPt;
105 }
106 else
107 {
109 << "Point perturbation crosses a surface. Not inserting."
110 << endl;
111 }
112 }
113
114 initialPoints.append(toPoint(midPoint));
115 }
116}
117
118
119// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
120
122(
123 const dictionary& initialPointsDict,
124 const Time& runTime,
125 Random& rndGen,
126 const conformationSurfaces& geometryToConformTo,
127 const cellShapeControl& cellShapeControls,
128 const autoPtr<backgroundMeshDecomposition>& decomposition
129)
130:
131 initialPointsMethod
132 (
133 typeName,
134 initialPointsDict,
135 runTime,
136 rndGen,
137 geometryToConformTo,
138 cellShapeControls,
139 decomposition
140 ),
141 randomiseInitialGrid_(detailsDict().get<Switch>("randomiseInitialGrid")),
142 randomPerturbationCoeff_
143 (
144 detailsDict().get<scalar>("randomPerturbationCoeff")
145 )
146{}
147
148
149// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
150
152{
153 // Loop over surface faces
154 const searchableSurfaces& surfaces = geometryToConformTo().geometry();
155 const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
156
157 const scalar maxRayLength(surfaces.bounds().mag());
158
159 // Initialise points list
160 label initialPointsSize = 0;
161 forAll(surfaces, surfI)
162 {
163 initialPointsSize += surfaces[surfI].size();
164 }
165
166 DynamicList<Vb::Point> initialPoints(initialPointsSize);
167
168 forAll(surfacesToConformTo, surfI)
169 {
170 const searchableSurface& s = surfaces[surfacesToConformTo[surfI]];
171
172 tmp<pointField> faceCentresTmp(s.coordinates());
173 const pointField& faceCentres = faceCentresTmp();
174
175 Info<< " Shoot rays from " << s.name() << nl
176 << " nRays = " << faceCentres.size() << endl;
177
178 forAll(faceCentres, fcI)
179 {
180 const Foam::point& fC = faceCentres[fcI];
181
182 if
183 (
185 && !decomposition().positionOnThisProcessor(fC)
186 )
187 {
188 continue;
189 }
190
191 const scalar pert
192 (
193 randomPerturbationCoeff_
194 *cellShapeControls().cellSize(fC)
195 );
196
197 pointIndexHit surfHitStart;
198 label hitSurfaceStart;
199
200 // Face centres should be on the surface so search distance can be
201 // small
203 (
204 fC,
205 sqr(pert),
206 surfHitStart,
207 hitSurfaceStart
208 );
209
210 vectorField normStart(1, vector::min);
212 (
213 hitSurfaceStart,
214 List<pointIndexHit>(1, surfHitStart),
215 normStart
216 );
217
218 pointIndexHit surfHitEnd;
219 label hitSurfaceEnd;
220
222 (
223 fC - normStart[0]*pert,
224 fC - normStart[0]*maxRayLength,
225 surfHitEnd,
226 hitSurfaceEnd
227 );
228
229 if (surfHitEnd.hit())
230 {
231 vectorField normEnd(1, vector::min);
233 (
234 hitSurfaceEnd,
235 List<pointIndexHit>(1, surfHitEnd),
236 normEnd
237 );
238
239 if ((normStart[0] & normEnd[0]) < 0)
240 {
241 line<point, point> l(fC, surfHitEnd.hitPoint());
242
243 if (Pstream::parRun())
244 {
245 // Clip the line in parallel
246 pointIndexHit procIntersection =
248 (
249 l.start(),
250 l.end()
251 );
252
253 if (procIntersection.hit())
254 {
255 l =
256 line<point, point>
257 (
258 l.start(),
259 procIntersection.hitPoint()
260 );
261 }
262 }
263
264 splitLine
265 (
266 l,
267 pert,
269 );
270 }
271 }
272 }
273 }
274
275 return initialPoints.shrink();
276}
277
278
279// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
void getNormal(const label hitSurface, const List< pointIndexHit > &surfHit, vectorField &normal) const
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
const searchableSurfaces & geometry() const
Return reference to the searchableSurfaces object containing all.
const labelList & surfaces() const
Return the surface indices.
void findSurfaceNearestIntersection(const point &start, const point &end, pointIndexHit &surfHit, label &hitSurface) const
Finding the nearestIntersection of the surface to start.
scalar minimumSurfaceDistanceCoeffSqr_
Only allow the placement of initial points that are within the.
const conformationSurfaces & geometryToConformTo() const
const cellShapeControl & cellShapeControls() const
const backgroundMeshDecomposition & decomposition() const
static const complex min
complex (-VGREAT,-VGREAT)
Definition: complex.H:292
virtual List< Vb::Point > initialPoints() const
Return the initial points for the conformalVoronoiMesh.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
engineTime & runTime
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointFrompoint toPoint(const Foam::point &p)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Field< vector > vectorField
Specialisation of Field<T> for vector.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Random rndGen
Definition: createFields.H:23