FreeFOAM The Cross-Platform CFD Toolkit
shellSurfaces.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-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 
27 #include "shellSurfaces.H"
28 #include <OpenFOAM/boundBox.H>
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 
38 namespace Foam
39 {
40 
41 template<>
42 const char*
43 NamedEnum<shellSurfaces::refineMode, 3>::
44 names[] =
45 {
46  "inside",
47  "outside",
48  "distance"
49 };
50 
51 const NamedEnum<shellSurfaces::refineMode, 3> shellSurfaces::refineModeNames_;
52 
53 } // End namespace Foam
54 
55 
56 
57 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
58 
59 void Foam::shellSurfaces::setAndCheckLevels
60 (
61  const label shellI,
62  const List<Tuple2<scalar, label> >& distLevels
63 )
64 {
65  if (modes_[shellI] != DISTANCE && distLevels.size() != 1)
66  {
68  (
69  "shellSurfaces::shellSurfaces"
70  "(const searchableSurfaces&, const dictionary&)"
71  ) << "For refinement mode "
72  << refineModeNames_[modes_[shellI]]
73  << " specify only one distance+level."
74  << " (its distance gets discarded)"
75  << exit(FatalError);
76  }
77  // Extract information into separate distance and level
78  distances_[shellI].setSize(distLevels.size());
79  levels_[shellI].setSize(distLevels.size());
80 
81  forAll(distLevels, j)
82  {
83  distances_[shellI][j] = distLevels[j].first();
84  levels_[shellI][j] = distLevels[j].second();
85 
86  // Check in incremental order
87  if (j > 0)
88  {
89  if
90  (
91  (distances_[shellI][j] <= distances_[shellI][j-1])
92  || (levels_[shellI][j] > levels_[shellI][j-1])
93  )
94  {
96  (
97  "shellSurfaces::shellSurfaces"
98  "(const searchableSurfaces&, const dictionary&)"
99  ) << "For refinement mode "
100  << refineModeNames_[modes_[shellI]]
101  << " : Refinement should be specified in order"
102  << " of increasing distance"
103  << " (and decreasing refinement level)." << endl
104  << "Distance:" << distances_[shellI][j]
105  << " refinementLevel:" << levels_[shellI][j]
106  << exit(FatalError);
107  }
108  }
109  }
110 
111  const searchableSurface& shell = allGeometry_[shells_[shellI]];
112 
113  if (modes_[shellI] == DISTANCE)
114  {
115  Info<< "Refinement level according to distance to "
116  << shell.name() << endl;
117  forAll(levels_[shellI], j)
118  {
119  Info<< " level " << levels_[shellI][j]
120  << " for all cells within " << distances_[shellI][j]
121  << " meter." << endl;
122  }
123  }
124  else
125  {
126  if (!allGeometry_[shells_[shellI]].hasVolumeType())
127  {
129  (
130  "shellSurfaces::shellSurfaces"
131  "(const searchableSurfaces&"
132  ", const PtrList<dictionary>&)"
133  ) << "Shell " << shell.name()
134  << " does not support testing for "
135  << refineModeNames_[modes_[shellI]] << endl
136  << "Probably it is not closed."
137  << exit(FatalError);
138  }
139 
140  if (modes_[shellI] == INSIDE)
141  {
142  Info<< "Refinement level " << levels_[shellI][0]
143  << " for all cells inside " << shell.name() << endl;
144  }
145  else
146  {
147  Info<< "Refinement level " << levels_[shellI][0]
148  << " for all cells outside " << shell.name() << endl;
149  }
150  }
151 }
152 
153 
154 // Specifically orient triSurfaces using a calculated point outside.
155 // Done since quite often triSurfaces not of consistent orientation which
156 // is (currently) necessary for sideness calculation
157 void Foam::shellSurfaces::orient()
158 {
159  // Determine outside point.
160  boundBox overallBb = boundBox::invertedBox;
161 
162  bool hasSurface = false;
163 
164  forAll(shells_, shellI)
165  {
166  const searchableSurface& s = allGeometry_[shells_[shellI]];
167 
168  if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
169  {
170  const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
171 
172  if (shell.triSurface::size())
173  {
174  const pointField& points = shell.points();
175 
176  hasSurface = true;
177 
178  boundBox shellBb(points[0], points[0]);
179  // Assume surface is compact!
180  for (label i = 0; i < points.size(); i++)
181  {
182  const point& pt = points[i];
183  shellBb.min() = min(shellBb.min(), pt);
184  shellBb.max() = max(shellBb.max(), pt);
185  }
186 
187  overallBb.min() = min(overallBb.min(), shellBb.min());
188  overallBb.max() = max(overallBb.max(), shellBb.max());
189  }
190  }
191  }
192 
193  if (hasSurface)
194  {
195  const point outsidePt = overallBb.max() + overallBb.span();
196 
197  //Info<< "Using point " << outsidePt << " to orient shells" << endl;
198 
199  forAll(shells_, shellI)
200  {
201  const searchableSurface& s = allGeometry_[shells_[shellI]];
202 
203  if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
204  {
205  triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
206  (
207  refCast<const triSurfaceMesh>(s)
208  );
209 
210  // Flip surface so outsidePt is outside.
211  bool anyFlipped = orientedSurface::orient
212  (
213  shell,
214  outsidePt,
215  true
216  );
217 
218  if (anyFlipped)
219  {
220  // orientedSurface will have done a clearOut of the surface.
221  // we could do a clearout of the triSurfaceMeshes::trees()
222  // but these aren't affected by orientation
223  // (except for cached
224  // sideness which should not be set at this point.
225  // !!Should check!)
226 
227  Info<< "shellSurfaces : Flipped orientation of surface "
228  << s.name()
229  << " so point " << outsidePt << " is outside." << endl;
230  }
231  }
232  }
233  }
234 }
235 
236 
237 // Find maximum level of a shell.
238 void Foam::shellSurfaces::findHigherLevel
239 (
240  const pointField& pt,
241  const label shellI,
242  labelList& maxLevel
243 ) const
244 {
245  const labelList& levels = levels_[shellI];
246 
247  if (modes_[shellI] == DISTANCE)
248  {
249  // Distance mode.
250 
251  const scalarField& distances = distances_[shellI];
252 
253  // Collect all those points that have a current maxLevel less than
254  // (any of) the shell. Also collect the furthest distance allowable
255  // to any shell with a higher level.
256 
257  pointField candidates(pt.size());
258  labelList candidateMap(pt.size());
259  scalarField candidateDistSqr(pt.size());
260  label candidateI = 0;
261 
262  forAll(maxLevel, pointI)
263  {
264  forAllReverse(levels, levelI)
265  {
266  if (levels[levelI] > maxLevel[pointI])
267  {
268  candidates[candidateI] = pt[pointI];
269  candidateMap[candidateI] = pointI;
270  candidateDistSqr[candidateI] = sqr(distances[levelI]);
271  candidateI++;
272  break;
273  }
274  }
275  }
276  candidates.setSize(candidateI);
277  candidateMap.setSize(candidateI);
278  candidateDistSqr.setSize(candidateI);
279 
280  // Do the expensive nearest test only for the candidate points.
281  List<pointIndexHit> nearInfo;
282  allGeometry_[shells_[shellI]].findNearest
283  (
284  candidates,
285  candidateDistSqr,
286  nearInfo
287  );
288 
289  // Update maxLevel
290  forAll(nearInfo, candidateI)
291  {
292  if (nearInfo[candidateI].hit())
293  {
294  // Check which level it actually is in.
295  label minDistI = findLower
296  (
297  distances,
298  mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
299  );
300 
301  label pointI = candidateMap[candidateI];
302 
303  // pt is inbetween shell[minDistI] and shell[minDistI+1]
304  maxLevel[pointI] = levels[minDistI+1];
305  }
306  }
307  }
308  else
309  {
310  // Inside/outside mode
311 
312  // Collect all those points that have a current maxLevel less than the
313  // shell.
314 
315  pointField candidates(pt.size());
316  labelList candidateMap(pt.size());
317  label candidateI = 0;
318 
319  forAll(maxLevel, pointI)
320  {
321  if (levels[0] > maxLevel[pointI])
322  {
323  candidates[candidateI] = pt[pointI];
324  candidateMap[candidateI] = pointI;
325  candidateI++;
326  }
327  }
328  candidates.setSize(candidateI);
329  candidateMap.setSize(candidateI);
330 
331  // Do the expensive nearest test only for the candidate points.
332  List<searchableSurface::volumeType> volType;
333  allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
334 
335  forAll(volType, i)
336  {
337  label pointI = candidateMap[i];
338 
339  if
340  (
341  (
342  modes_[shellI] == INSIDE
343  && volType[i] == searchableSurface::INSIDE
344  )
345  || (
346  modes_[shellI] == OUTSIDE
347  && volType[i] == searchableSurface::OUTSIDE
348  )
349  )
350  {
351  maxLevel[pointI] = levels[0];
352  }
353  }
354  }
355 }
356 
357 
358 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
359 
361 (
362  const searchableSurfaces& allGeometry,
363  const PtrList<dictionary>& shellDicts
364 )
365 :
366  allGeometry_(allGeometry)
367 {
368  shells_.setSize(shellDicts.size());
369  modes_.setSize(shellDicts.size());
370  distances_.setSize(shellDicts.size());
371  levels_.setSize(shellDicts.size());
372 
373  forAll(shellDicts, shellI)
374  {
375  const dictionary& dict = shellDicts[shellI];
376  const word name = dict.lookup("name");
377  const word type = dict.lookup("type");
378 
379  shells_[shellI] = allGeometry_.findSurfaceID(name);
380 
381  if (shells_[shellI] == -1)
382  {
384  (
385  "shellSurfaces::shellSurfaces"
386  "(const searchableSurfaces&, const PtrList<dictionary>&)"
387  ) << "No surface called " << name << endl
388  << "Valid surfaces are " << allGeometry_.names()
389  << exit(FatalError);
390  }
391 
392  modes_[shellI] = refineModeNames_.read(dict.lookup("refineMode"));
393 
394  // Read pairs of distance+level
395  setAndCheckLevels(shellI, dict.lookup("levels"));
396  }
397 
398  // Orient shell surfaces before any searching is done. Note that this
399  // only needs to be done for inside or outside. Orienting surfaces
400  // constructs lots of addressing which we want to avoid.
401  orient();
402 }
403 
404 
406 (
407  const searchableSurfaces& allGeometry,
408  const dictionary& shellsDict
409 )
410 :
411  allGeometry_(allGeometry)
412 {
413  shells_.setSize(shellsDict.size());
414  modes_.setSize(shellsDict.size());
415  distances_.setSize(shellsDict.size());
416  levels_.setSize(shellsDict.size());
417 
418  label shellI = 0;
419  forAllConstIter(dictionary, shellsDict, iter)
420  {
421  shells_[shellI] = allGeometry_.findSurfaceID(iter().keyword());
422 
423  if (shells_[shellI] == -1)
424  {
426  (
427  "shellSurfaces::shellSurfaces"
428  "(const searchableSurfaces&, const dictionary>&"
429  ) << "No surface called " << iter().keyword() << endl
430  << "Valid surfaces are " << allGeometry_.names()
431  << exit(FatalError);
432  }
433  const dictionary& dict = shellsDict.subDict(iter().keyword());
434 
435  modes_[shellI] = refineModeNames_.read(dict.lookup("mode"));
436 
437  // Read pairs of distance+level
438  setAndCheckLevels(shellI, dict.lookup("levels"));
439 
440  shellI++;
441  }
442 
443  // Orient shell surfaces before any searching is done. Note that this
444  // only needs to be done for inside or outside. Orienting surfaces
445  // constructs lots of addressing which we want to avoid.
446  orient();
447 }
448 
449 
450 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
451 
452 // Highest shell level
453 Foam::label Foam::shellSurfaces::maxLevel() const
454 {
455  label overallMax = 0;
456  forAll(levels_, shellI)
457  {
458  overallMax = max(overallMax, max(levels_[shellI]));
459  }
460  return overallMax;
461 }
462 
463 
464 void Foam::shellSurfaces::findHigherLevel
465 (
466  const pointField& pt,
467  const labelList& ptLevel,
468  labelList& maxLevel
469 ) const
470 {
471  // Maximum level of any shell. Start off with level of point.
472  maxLevel = ptLevel;
473 
474  forAll(shells_, shellI)
475  {
476  findHigherLevel(pt, shellI, maxLevel);
477  }
478 }
479 
480 
481 // ************************ vim: set sw=4 sts=4 et: ************************ //