FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2011 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 
28 #include <OpenFOAM/Time.H>
29 #include <OpenFOAM/ListOps.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
37 addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
38 
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
45 (
46  const point& start,
47  const point& end
48 ) const
49 {
50  Pair<vector> offsets(vector::zero, vector::zero);
51 
52  vector n(end-start);
53 
54  scalar magN = mag(n);
55 
56  if (magN > SMALL)
57  {
58  n /= magN;
59 
60  // Do first offset vector. Is the coordinate axes with the smallest
61  // component along the vector n.
62  scalar minMag = GREAT;
63  direction minCmpt = 0;
64 
65  for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
66  {
67  if (mag(n[cmpt]) < minMag)
68  {
69  minMag = mag(n[cmpt]);
70  minCmpt = cmpt;
71  }
72  }
73 
74  offsets[0][minCmpt] = 1.0;
75  // Orthonormalise
76  offsets[0] -= n[minCmpt]*n;
77  offsets[0] /= mag(offsets[0]);
78  // Do second offset vector perp to original edge and first offset vector
79  offsets[1] = n ^ offsets[0];
80 
81  // Scale
82  offsets[0] *= gap_;
83  offsets[1] *= gap_;
84  }
85 
86  return offsets;
87 }
88 
89 
90 void Foam::searchableSurfaceWithGaps::offsetVecs
91 (
92  const pointField& start,
93  const pointField& end,
94  pointField& offset0,
95  pointField& offset1
96 ) const
97 {
98  offset0.setSize(start.size());
99  offset1.setSize(start.size());
100 
101  forAll(start, i)
102  {
103  const Pair<vector> offsets(offsetVecs(start[i], end[i]));
104  offset0[i] = offsets[0];
105  offset1[i] = offsets[1];
106  }
107 }
108 
109 
110 Foam::label Foam::searchableSurfaceWithGaps::countMisses
111 (
112  const List<pointIndexHit>& info,
113  labelList& missMap
114 )
115 {
116  label nMiss = 0;
117  forAll(info, i)
118  {
119  if (!info[i].hit())
120  {
121  nMiss++;
122  }
123  }
124 
125  missMap.setSize(nMiss);
126  nMiss = 0;
127 
128  forAll(info, i)
129  {
130  if (!info[i].hit())
131  {
132  missMap[nMiss++] = i;
133  }
134  }
135 
136  return nMiss;
137 }
138 
139 
140 // Anything not a hit in both counts as a hit
141 Foam::label Foam::searchableSurfaceWithGaps::countMisses
142 (
143  const List<pointIndexHit>& plusInfo,
144  const List<pointIndexHit>& minInfo,
145  labelList& missMap
146 )
147 {
148  label nMiss = 0;
149  forAll(plusInfo, i)
150  {
151  if (!plusInfo[i].hit() || !minInfo[i].hit())
152  {
153  nMiss++;
154  }
155  }
156 
157  missMap.setSize(nMiss);
158  nMiss = 0;
159 
160  forAll(plusInfo, i)
161  {
162  if (!plusInfo[i].hit() || !minInfo[i].hit())
163  {
164  missMap[nMiss++] = i;
165  }
166  }
167 
168  return nMiss;
169 }
170 
171 
172 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 
174 Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
175 (
176  const IOobject& io,
177  const dictionary& dict
178 )
179 :
180  searchableSurface(io),
181  gap_(readScalar(dict.lookup("gap"))),
182  subGeom_(1)
183 {
184  const word subGeomName(dict.lookup("surface"));
185 
186  const searchableSurface& s =
187  io.db().lookupObject<searchableSurface>(subGeomName);
188 
189  subGeom_.set(0, &const_cast<searchableSurface&>(s));
190 }
191 
192 
193 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
194 
196 {}
197 
198 
199 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
200 
202 (
203  const pointField& start,
204  const pointField& end,
205  List<pointIndexHit>& info
206 ) const
207 {
208 
209  // Test with unperturbed vectors
210  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211 
212  surface().findLine(start, end, info);
213 
214  // Count number of misses. Determine map
215  labelList compactMap;
216  label nMiss = countMisses(info, compactMap);
217 
218  if (returnReduce(nMiss, sumOp<label>()) > 0)
219  {
220  //Pout<< "** retesting with offset0 " << nMiss << " misses out of "
221  // << start.size() << endl;
222 
223  // extract segments according to map
224  pointField compactStart(start, compactMap);
225  pointField compactEnd(end, compactMap);
226 
227  // Calculate offset vector
228  pointField offset0, offset1;
229  offsetVecs
230  (
231  compactStart,
232  compactEnd,
233  offset0,
234  offset1
235  );
236 
237  // Test with offset0 perturbed vectors
238  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
239 
240  // test in pairs: only if both perturbations hit something
241  // do we accept the hit.
242 
243  const vectorField smallVec(SMALL*(compactEnd-compactStart));
244 
245  List<pointIndexHit> plusInfo;
246  surface().findLine
247  (
248  compactStart+offset0-smallVec,
249  compactEnd+offset0+smallVec,
250  plusInfo
251  );
252  List<pointIndexHit> minInfo;
253  surface().findLine
254  (
255  compactStart-offset0-smallVec,
256  compactEnd-offset0+smallVec,
257  minInfo
258  );
259 
260  // Extract any hits
261  forAll(plusInfo, i)
262  {
263  if (plusInfo[i].hit() && minInfo[i].hit())
264  {
265  info[compactMap[i]] = plusInfo[i];
266  info[compactMap[i]].rawPoint() -= offset0[i];
267  }
268  }
269 
270  labelList plusMissMap;
271  nMiss = countMisses(plusInfo, minInfo, plusMissMap);
272 
273  if (returnReduce(nMiss, sumOp<label>()) > 0)
274  {
275  //Pout<< "** retesting with offset1 " << nMiss << " misses out of "
276  // << start.size() << endl;
277 
278  // Test with offset1 perturbed vectors
279  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
280 
281  // Extract (inplace possible because of order)
282  forAll(plusMissMap, i)
283  {
284  label mapI = plusMissMap[i];
285  compactStart[i] = compactStart[mapI];
286  compactEnd[i] = compactEnd[mapI];
287  compactMap[i] = compactMap[mapI];
288  offset0[i] = offset0[mapI];
289  offset1[i] = offset1[mapI];
290  }
291  compactStart.setSize(plusMissMap.size());
292  compactEnd.setSize(plusMissMap.size());
293  compactMap.setSize(plusMissMap.size());
294  offset0.setSize(plusMissMap.size());
295  offset1.setSize(plusMissMap.size());
296 
297  const vectorField smallVec(SMALL*(compactEnd-compactStart));
298 
299  surface().findLine
300  (
301  compactStart+offset1-smallVec,
302  compactEnd+offset1+smallVec,
303  plusInfo
304  );
305  surface().findLine
306  (
307  compactStart-offset1-smallVec,
308  compactEnd-offset1+smallVec,
309  minInfo
310  );
311 
312  // Extract any hits
313  forAll(plusInfo, i)
314  {
315  if (plusInfo[i].hit() && minInfo[i].hit())
316  {
317  info[compactMap[i]] = plusInfo[i];
318  info[compactMap[i]].rawPoint() -= offset1[i];
319  }
320  }
321  }
322  }
323 }
324 
325 
327 (
328  const pointField& start,
329  const pointField& end,
330  List<pointIndexHit>& info
331 ) const
332 {
333  // To be done ...
334  findLine(start, end, info);
335 }
336 
337 
339 (
340  const pointField& start,
341  const pointField& end,
342  List<List<pointIndexHit> >& info
343 ) const
344 {
345  // To be done. Assume for now only one intersection.
346  List<pointIndexHit> nearestInfo;
347  findLine(start, end, nearestInfo);
348 
349  info.setSize(start.size());
350  forAll(info, pointI)
351  {
352  if (nearestInfo[pointI].hit())
353  {
354  info[pointI].setSize(1);
355  info[pointI][0] = nearestInfo[pointI];
356  }
357  else
358  {
359  info[pointI].clear();
360  }
361  }
362 }
363 
364 
365 // ************************ vim: set sw=4 sts=4 et: ************************ //