searchableSurfaceWithGaps.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) 2011-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
30#include "Time.H"
31#include "ListOps.H"
32
33// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34
35namespace Foam
36{
39 (
42 dict
43 );
44}
45
46
47// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48
49Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
50(
51 const point& start,
52 const point& end
53) const
54{
55 Pair<vector> offsets(Zero, Zero);
56
57 vector n(end-start);
58
59 scalar magN = mag(n);
60
61 if (magN > SMALL)
62 {
63 n /= magN;
64
65 // Do first offset vector. Is the coordinate axes with the smallest
66 // component along the vector n.
67 scalar minMag = GREAT;
68 direction minCmpt = 0;
69
70 for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
71 {
72 if (mag(n[cmpt]) < minMag)
73 {
74 minMag = mag(n[cmpt]);
75 minCmpt = cmpt;
76 }
77 }
78
79 offsets[0][minCmpt] = 1.0;
80 // Orthonormalise
81 offsets[0] -= n[minCmpt]*n;
82 offsets[0] /= mag(offsets[0]);
83 // Do second offset vector perp to original edge and first offset vector
84 offsets[1] = n ^ offsets[0];
85
86 // Scale
87 offsets[0] *= gap_;
88 offsets[1] *= gap_;
89 }
90
91 return offsets;
92}
93
94
95void Foam::searchableSurfaceWithGaps::offsetVecs
96(
97 const pointField& start,
98 const pointField& end,
99 pointField& offset0,
100 pointField& offset1
101) const
102{
103 offset0.setSize(start.size());
104 offset1.setSize(start.size());
105
106 forAll(start, i)
107 {
108 const Pair<vector> offsets(offsetVecs(start[i], end[i]));
109 offset0[i] = offsets[0];
110 offset1[i] = offsets[1];
111 }
112}
113
114
115Foam::label Foam::searchableSurfaceWithGaps::countMisses
116(
117 const List<pointIndexHit>& info,
118 labelList& missMap
119)
120{
121 label nMiss = 0;
122 forAll(info, i)
123 {
124 if (!info[i].hit())
125 {
126 nMiss++;
127 }
128 }
129
130 missMap.setSize(nMiss);
131 nMiss = 0;
132
133 forAll(info, i)
134 {
135 if (!info[i].hit())
136 {
137 missMap[nMiss++] = i;
138 }
139 }
140
141 return nMiss;
142}
143
144
145// Anything not a hit in both counts as a hit
146Foam::label Foam::searchableSurfaceWithGaps::countMisses
147(
148 const List<pointIndexHit>& plusInfo,
149 const List<pointIndexHit>& minInfo,
150 labelList& missMap
151)
152{
153 label nMiss = 0;
154 forAll(plusInfo, i)
155 {
156 if (!plusInfo[i].hit() || !minInfo[i].hit())
157 {
158 nMiss++;
159 }
160 }
161
162 missMap.setSize(nMiss);
163 nMiss = 0;
164
165 forAll(plusInfo, i)
166 {
167 if (!plusInfo[i].hit() || !minInfo[i].hit())
168 {
169 missMap[nMiss++] = i;
170 }
171 }
172
173 return nMiss;
174}
175
176
177// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
178
180(
181 const IOobject& io,
182 const dictionary& dict
183)
184:
186 gap_(dict.get<scalar>("gap")),
187 subGeom_(1)
188{
189 const word subGeomName(dict.get<word>("surface"));
190
191 subGeom_.set
192 (
193 0,
194 io.db().getObjectPtr<searchableSurface>(subGeomName)
195 );
196
197 bounds() = subGeom_[0].bounds();
198}
199
200
201// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
202
204(
205 const pointField& start,
206 const pointField& end,
208) const
209{
210
211 // Test with unperturbed vectors
212 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213
214 surface().findLine(start, end, info);
215
216 // Count number of misses. Determine map
217 labelList compactMap;
218 label nMiss = countMisses(info, compactMap);
219
220 if (returnReduce(nMiss, sumOp<label>()) > 0)
221 {
222 //Pout<< "** retesting with offset0 " << nMiss << " misses out of "
223 // << start.size() << endl;
224
225 // extract segments according to map
226 pointField compactStart(start, compactMap);
227 pointField compactEnd(end, compactMap);
228
229 // Calculate offset vector
230 pointField offset0, offset1;
231 offsetVecs
232 (
233 compactStart,
234 compactEnd,
235 offset0,
236 offset1
237 );
238
239 // Test with offset0 perturbed vectors
240 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
241
242 // test in pairs: only if both perturbations hit something
243 // do we accept the hit.
244
245 const vectorField smallVec(1e-6*(compactEnd-compactStart));
246
247 List<pointIndexHit> plusInfo;
248 surface().findLine
249 (
250 compactStart+offset0-smallVec,
251 compactEnd+offset0+smallVec,
252 plusInfo
253 );
254 List<pointIndexHit> minInfo;
255 surface().findLine
256 (
257 compactStart-offset0-smallVec,
258 compactEnd-offset0+smallVec,
259 minInfo
260 );
261
262 // Extract any hits
263 forAll(plusInfo, i)
264 {
265 if (plusInfo[i].hit() && minInfo[i].hit())
266 {
267 info[compactMap[i]] = plusInfo[i];
268 info[compactMap[i]].rawPoint() -= offset0[i];
269 }
270 }
271
272 labelList plusMissMap;
273 nMiss = countMisses(plusInfo, minInfo, plusMissMap);
274
275 if (returnReduce(nMiss, sumOp<label>()) > 0)
276 {
277 //Pout<< "** retesting with offset1 " << nMiss << " misses out of "
278 // << start.size() << endl;
279
280 // Test with offset1 perturbed vectors
281 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282
283 // Extract (inplace possible because of order)
284 forAll(plusMissMap, i)
285 {
286 label mapI = plusMissMap[i];
287 compactStart[i] = compactStart[mapI];
288 compactEnd[i] = compactEnd[mapI];
289 compactMap[i] = compactMap[mapI];
290 offset0[i] = offset0[mapI];
291 offset1[i] = offset1[mapI];
292 }
293 compactStart.setSize(plusMissMap.size());
294 compactEnd.setSize(plusMissMap.size());
295 compactMap.setSize(plusMissMap.size());
296 offset0.setSize(plusMissMap.size());
297 offset1.setSize(plusMissMap.size());
298
299 const vectorField smallVec(1e-6*(compactEnd-compactStart));
300
301 surface().findLine
302 (
303 compactStart+offset1-smallVec,
304 compactEnd+offset1+smallVec,
305 plusInfo
306 );
307 surface().findLine
308 (
309 compactStart-offset1-smallVec,
310 compactEnd-offset1+smallVec,
311 minInfo
312 );
313
314 // Extract any hits
315 forAll(plusInfo, i)
316 {
317 if (plusInfo[i].hit() && minInfo[i].hit())
318 {
319 info[compactMap[i]] = plusInfo[i];
320 info[compactMap[i]].rawPoint() -= offset1[i];
321 }
322 }
323 }
324 }
325}
326
327
329(
330 const pointField& start,
331 const pointField& end,
333) const
334{
335 // To be done ...
336 findLine(start, end, info);
337}
338
339
341(
342 const pointField& start,
343 const pointField& end,
345) const
346{
347 // To be done. Assume for now only one intersection.
348 List<pointIndexHit> nearestInfo;
349 findLine(start, end, nearestInfo);
350
351 info.setSize(start.size());
352 forAll(info, pointi)
353 {
354 if (nearestInfo[pointi].hit())
355 {
356 info[pointi].setSize(1);
357 info[pointi][0] = nearestInfo[pointi];
358 }
359 else
360 {
361 info[pointi].clear();
362 }
363 }
364}
365
366
367// ************************************************************************* //
Various functions to operate on Lists.
label n
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const objectRegistry & db() const noexcept
Return the local objectRegistry.
Definition: IOobject.C:500
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
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
const T * set(const label i) const
Definition: UPtrList.H:248
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Type * getObjectPtr(const word &name, const bool recursive=false) const
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
searchableSurface using multiple slightly shifted underlying surfaces to make sure pierces don't go t...
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit > > &) const
Get all intersections in order from start to end.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
virtual const boundBox & bounds() const
Return const reference to boundBox.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333