All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
shellSurfaces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 "shellSurfaces.H"
27 #include "refinementSurfaces.H"
28 #include "searchableSurfaces.H"
29 #include "orientedSurface.H"
30 #include "volumeType.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 template<>
38 const char*
39 NamedEnum<shellSurfaces::refineMode, 5>::
40 names[] =
41 {
42  "inside",
43  "outside",
44  "distance",
45  "insideSpan",
46  "outsideSpan"
47 };
48 
49 }
50 
52  Foam::shellSurfaces::refineModeNames_;
53 
54 
55 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 
57 void Foam::shellSurfaces::setAndCheckLevels
58 (
59  const label shelli,
60  const List<Tuple2<scalar, label>>& distLevels
61 )
62 {
63  if (modes_[shelli] != refineMode::distance && distLevels.size() != 1)
64  {
66  << "For refinement mode "
67  << refineModeNames_[modes_[shelli]]
68  << " specify only one distance+level."
69  << " (its distance gets discarded)"
70  << exit(FatalError);
71  }
72 
73  // Extract information into separate distance and level
74  distances_[shelli].setSize(distLevels.size());
75  levels_[shelli].setSize(distLevels.size());
76 
77  forAll(distLevels, j)
78  {
79  distances_[shelli][j] = distLevels[j].first();
80  levels_[shelli][j] = distLevels[j].second();
81 
82  // Check in incremental order
83  if (j > 0)
84  {
85  if
86  (
87  (distances_[shelli][j] <= distances_[shelli][j-1])
88  || (levels_[shelli][j] > levels_[shelli][j-1])
89  )
90  {
92  << "For refinement mode "
93  << refineModeNames_[modes_[shelli]]
94  << " : Refinement should be specified in order"
95  << " of increasing distance"
96  << " (and decreasing refinement level)." << endl
97  << "Distance:" << distances_[shelli][j]
98  << " refinementLevel:" << levels_[shelli][j]
99  << exit(FatalError);
100  }
101  }
102  }
103 
104  const searchableSurface& shell = allGeometry_[shells_[shelli]];
105 
106  if (modes_[shelli] == refineMode::distance)
107  {
108  Info<< "Refinement level according to distance to "
109  << shell.name() << endl;
110 
111  forAll(levels_[shelli], j)
112  {
113  Info<< " level " << levels_[shelli][j]
114  << " for all cells within " << distances_[shelli][j]
115  << " metre." << endl;
116  }
117  }
118  else
119  {
120  if (!allGeometry_[shells_[shelli]].hasVolumeType())
121  {
123  << "Shell " << shell.name()
124  << " does not support testing for "
125  << refineModeNames_[modes_[shelli]] << endl
126  << "Probably it is not closed."
127  << exit(FatalError);
128  }
129 
130  if (modes_[shelli] == refineMode::inside)
131  {
132  Info<< "Refinement level " << levels_[shelli][0]
133  << " for all cells inside " << shell.name() << endl;
134  }
135  else
136  {
137  Info<< "Refinement level " << levels_[shelli][0]
138  << " for all cells outside " << shell.name() << endl;
139  }
140  }
141 }
142 
143 
144 void Foam::shellSurfaces::orient()
145 {
146  // Determine outside point.
147  boundBox overallBb = boundBox::invertedBox;
148 
149  bool hasSurface = false;
150 
151  forAll(shells_, shelli)
152  {
153  const searchableSurface& s = allGeometry_[shells_[shelli]];
154 
155  if (modes_[shelli] != refineMode::distance && isA<triSurfaceMesh>(s))
156  {
157  const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
158 
159  if (shell.triSurface::size())
160  {
161  const pointField& points = shell.points();
162 
163  hasSurface = true;
164 
165  boundBox shellBb(points[0], points[0]);
166  // Assume surface is compact!
167  forAll(points, i)
168  {
169  const point& pt = points[i];
170  shellBb.min() = min(shellBb.min(), pt);
171  shellBb.max() = max(shellBb.max(), pt);
172  }
173 
174  overallBb.min() = min(overallBb.min(), shellBb.min());
175  overallBb.max() = max(overallBb.max(), shellBb.max());
176  }
177  }
178  }
179 
180  if (hasSurface)
181  {
182  const point outsidePt = overallBb.max() + overallBb.span();
183 
184  // Info<< "Using point " << outsidePt << " to orient shells" << endl;
185 
186  forAll(shells_, shelli)
187  {
188  const searchableSurface& s = allGeometry_[shells_[shelli]];
189 
190  if
191  (
192  modes_[shelli] != refineMode::distance
193  && isA<triSurfaceMesh>(s)
194  )
195  {
196  triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
197  (
198  refCast<const triSurfaceMesh>(s)
199  );
200 
201  // Flip surface so outsidePt is outside.
202  bool anyFlipped = orientedSurface::orient
203  (
204  shell,
205  outsidePt,
206  true
207  );
208 
209  if (anyFlipped)
210  {
211  // orientedSurface will have done a clearOut of the surface.
212  // we could do a clearout of the triSurfaceMeshes::trees()
213  // but these aren't affected by orientation
214  // (except for cached
215  // sideness which should not be set at this point.
216  // !!Should check!)
217 
218  Info<< "shellSurfaces : Flipped orientation of surface "
219  << s.name()
220  << " so point " << outsidePt << " is outside." << endl;
221  }
222  }
223  }
224  }
225 }
226 
227 
228 Foam::scalar Foam::shellSurfaces::interpolate
229 (
230  const triSurfaceMesh& tsm,
231  const triSurfacePointScalarField& closeness,
232  const point& pt,
233  const label index
234 ) const
235 {
236  const barycentric2D bary
237  (
239  (
240  tsm.points(),
241  tsm.triSurface::operator[](index)
242  ).pointToBarycentric(pt)
243  );
244 
245  const labelledTri& lf = tsm.localFaces()[index];
246 
247  return closeness[lf[0]]*bary[0]
248  + closeness[lf[1]]*bary[1]
249  + closeness[lf[2]]*bary[2];
250 }
251 
252 
253 void Foam::shellSurfaces::findHigherLevel
254 (
255  const pointField& pt,
256  const label shelli,
257  const scalar level0EdgeLength,
258  labelList& maxLevel
259 ) const
260 {
261  const labelList& levels = levels_[shelli];
262 
263  if (modes_[shelli] == refineMode::distance)
264  {
265  // Distance mode.
266 
267  const scalarField& distances = distances_[shelli];
268 
269  // Collect all those points that have a current maxLevel less than
270  // (any of) the shell. Also collect the furthest distance allowable
271  // to any shell with a higher level.
272 
273  pointField candidates(pt.size());
274  labelList candidateMap(pt.size());
275  scalarField candidateDistSqr(pt.size());
276  label candidateI = 0;
277 
278  forAll(maxLevel, pointi)
279  {
280  forAllReverse(levels, levelI)
281  {
282  if (levels[levelI] > maxLevel[pointi])
283  {
284  candidates[candidateI] = pt[pointi];
285  candidateMap[candidateI] = pointi;
286  candidateDistSqr[candidateI] = sqr(distances[levelI]);
287  candidateI++;
288  break;
289  }
290  }
291  }
292  candidates.setSize(candidateI);
293  candidateMap.setSize(candidateI);
294  candidateDistSqr.setSize(candidateI);
295 
296  // Do the expensive nearest test only for the candidate points.
297  List<pointIndexHit> nearInfo;
298  allGeometry_[shells_[shelli]].findNearest
299  (
300  candidates,
301  candidateDistSqr,
302  nearInfo
303  );
304 
305  // Update maxLevel
306  forAll(nearInfo, candidateI)
307  {
308  if (nearInfo[candidateI].hit())
309  {
310  // Check which level it actually is in.
311  label minDistI = findLower
312  (
313  distances,
314  mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
315  );
316 
317  label pointi = candidateMap[candidateI];
318 
319  // pt is in between shell[minDistI] and shell[minDistI+1]
320  maxLevel[pointi] = levels[minDistI+1];
321  }
322  }
323  }
324  else if
325  (
326  modes_[shelli] == refineMode::insideSpan
327  || modes_[shelli] == refineMode::outsideSpan
328  )
329  {
330  const triSurfaceMesh& tsm =
331  refCast<const triSurfaceMesh>(allGeometry_[shells_[shelli]]);
332 
333  // Collect all those points that have a current maxLevel less than
334  // the maximum and furthest distance allowable for the shell.
335 
336  pointField candidates(pt.size());
337  labelList candidateMap(pt.size());
338  scalarField candidateDistSqr(pt.size());
339  label candidateI = 0;
340 
341  forAll(pt, pointi)
342  {
343  if (levels[0] > maxLevel[pointi])
344  {
345  candidates[candidateI] = pt[pointi];
346  candidateMap[candidateI] = pointi;
347  candidateDistSqr[candidateI] = sqr(distances_[shelli][0]);
348  candidateI++;
349  }
350  }
351  candidates.setSize(candidateI);
352  candidateMap.setSize(candidateI);
353  candidateDistSqr.setSize(candidateI);
354 
355  // Do the expensive nearest test only for the candidate points.
356  List<pointIndexHit> nearInfo;
357  tsm.findNearest
358  (
359  candidates,
360  candidateDistSqr,
361  nearInfo
362  );
363 
364  // Minimum span for the maximum level specified
365  const scalar minSpan
366  (
367  cellsAcrossSpan_[shelli]*level0EdgeLength/(1 << (levels[0] - 1))
368  );
369 
370  // Update maxLevel
371  forAll(nearInfo, candidateI)
372  {
373  if (nearInfo[candidateI].hit())
374  {
375  const scalar span
376  (
378  (
379  tsm,
380  closeness_[shelli],
381  nearInfo[candidateI].rawPoint(),
382  nearInfo[candidateI].index()
383  )
384  );
385 
386  if (span > minSpan)
387  {
388  const label level
389  (
390  log2(cellsAcrossSpan_[shelli]*level0EdgeLength/span) + 1
391  );
392 
393  maxLevel[candidateMap[candidateI]] = min(levels[0], level);
394  }
395  else
396  {
397  maxLevel[candidateMap[candidateI]] = levels[0];
398  }
399  }
400  }
401  }
402  else
403  {
404  // Inside/outside mode
405 
406  // Collect all those points that have a current maxLevel less than the
407  // shell.
408 
409  pointField candidates(pt.size());
410  labelList candidateMap(pt.size());
411  label candidateI = 0;
412 
413  forAll(maxLevel, pointi)
414  {
415  if (levels[0] > maxLevel[pointi])
416  {
417  candidates[candidateI] = pt[pointi];
418  candidateMap[candidateI] = pointi;
419  candidateI++;
420  }
421  }
422  candidates.setSize(candidateI);
423  candidateMap.setSize(candidateI);
424 
425  // Do the expensive nearest test only for the candidate points.
426  List<volumeType> volType;
427  allGeometry_[shells_[shelli]].getVolumeType(candidates, volType);
428 
429  forAll(volType, i)
430  {
431  label pointi = candidateMap[i];
432 
433  if
434  (
435  (
436  modes_[shelli] == refineMode::inside
437  && volType[i] == volumeType::inside
438  )
439  || (
440  modes_[shelli] == refineMode::outside
441  && volType[i] == volumeType::outside
442  )
443  )
444  {
445  maxLevel[pointi] = levels[0];
446  }
447  }
448  }
449 }
450 
451 
452 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
453 
455 (
456  const searchableSurfaces& allGeometry,
457  const dictionary& shellsDict
458 )
459 :
460  allGeometry_(allGeometry)
461 {
462  // Wildcard specification : loop over all surfaces and try to find a match.
463 
464  // Count number of shells.
465  label shelli = 0;
466  forAll(allGeometry.names(), geomI)
467  {
468  const word& geomName = allGeometry_.names()[geomI];
469 
470  if (shellsDict.found(geomName))
471  {
472  shelli++;
473  }
474  }
475 
476  // Size lists
477  shells_.setSize(shelli);
478  modes_.setSize(shelli);
479  distances_.setSize(shelli);
480  levels_.setSize(shelli);
481  cellsAcrossSpan_.setSize(shelli);
482  closeness_.setSize(shelli);
483 
484  HashSet<word> unmatchedKeys(shellsDict.toc());
485  shelli = 0;
486 
487  forAll(allGeometry_.names(), geomI)
488  {
489  const word& geomName = allGeometry_.names()[geomI];
490 
491  const entry* ePtr = shellsDict.lookupEntryPtr(geomName, false, true);
492 
493  if (ePtr)
494  {
495  const dictionary& dict = ePtr->dict();
496  unmatchedKeys.erase(ePtr->keyword());
497 
498  shells_[shelli] = geomI;
499  modes_[shelli] = refineModeNames_.read(dict.lookup("mode"));
500 
501  // Read pairs of distance+level
502  setAndCheckLevels(shelli, dict.lookup("levels"));
503 
504  if
505  (
506  modes_[shelli] == refineMode::insideSpan
507  || modes_[shelli] == refineMode::outsideSpan
508  )
509  {
510  const searchableSurface& surface = allGeometry_[geomI];
511 
512  if (isA<triSurfaceMesh>(surface))
513  {
514  dict.lookup("cellsAcrossSpan") >> cellsAcrossSpan_[shelli];
515 
516  const triSurfaceMesh& tsm =
517  refCast<const triSurfaceMesh>(surface);
518 
519  closeness_.set
520  (
521  shelli,
523  (
524  IOobject
525  (
526  surface.name()
527  + ".closeness."
528  + (
529  modes_[shelli] == refineMode::insideSpan
530  ? "internal"
531  : "external"
532  )
533  + "PointCloseness",
534  surface.searchableSurface::time().constant(),
535  "triSurface",
536  surface.searchableSurface::time(),
538  ),
539  tsm,
540  dimLength,
541  true
542  )
543  );
544  }
545  else
546  {
547  FatalIOErrorInFunction(shellsDict)
548  << "Surface " << surface.name()
549  << " is not a triSurface as required by"
550  " refinement modes "
551  << refineModeNames_[refineMode::insideSpan]
552  << " and " << refineModeNames_[refineMode::outsideSpan]
553  << exit(FatalIOError);
554  }
555  }
556 
557  shelli++;
558  }
559  }
560 
561  if (unmatchedKeys.size() > 0)
562  {
564  (
565  shellsDict
566  ) << "Not all entries in refinementRegions dictionary were used."
567  << " The following entries were not used : "
568  << unmatchedKeys.sortedToc()
569  << endl;
570  }
571 
572  // Orient shell surfaces before any searching is done. Note that this
573  // only needs to be done for inside or outside. Orienting surfaces
574  // constructs lots of addressing which we want to avoid.
575  orient();
576 }
577 
578 
579 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
580 
582 {
583  label overallMax = 0;
584 
585  forAll(levels_, shelli)
586  {
587  overallMax = max(overallMax, max(levels_[shelli]));
588  }
589 
590  return overallMax;
591 }
592 
593 
594 void Foam::shellSurfaces::findHigherLevel
595 (
596  const pointField& pt,
597  const labelList& ptLevel,
598  const scalar level0EdgeLength,
599  labelList& maxLevel
600 ) const
601 {
602  // Maximum level of any shell. Start off with level of point.
603  maxLevel = ptLevel;
604 
605  forAll(shells_, shelli)
606  {
607  findHigherLevel(pt, shelli, level0EdgeLength, maxLevel);
608  }
609 }
610 
611 
612 // ************************************************************************* //
shellSurfaces(const searchableSurfaces &allGeometry, const dictionary &shellsDict)
Construct from geometry and dictionary.
A HashTable with keys but without contents.
Definition: HashSet.H:59
dictionary dict
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:667
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
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
const word & name() const
Return name.
Definition: IOobject.H:303
const entry * lookupEntryPtr(const word &, bool recursive, bool patternMatch) const
Find and return an entry data stream pointer if present.
Definition: dictionary.C:706
const keyType & keyword() const
Return keyword.
Definition: entry.H:123
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
barycentric2D pointToBarycentric(const point &pt) const
Calculate the barycentric coordinates from the given point.
Definition: triangleI.H:259
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label log2(label i)
Return the log base 2 by successive bit-shifting of the given label.
Definition: label.H:79
A 2-tuple for storing two objects of different types.
Definition: HashTable.H:65
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
virtual tmp< pointField > points() const
Get the points that define the surface.
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool orient(triSurface &, const point &, const bool orientOutside)
Flip faces such that normals are consistent with point:
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
virtual const dictionary & dict() const =0
Return dictionary if this entry is a dictionary.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
label maxLevel() const
Highest shell level.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
bool erase(T *p)
Remove the specified element from the list and delete it.
Definition: ILList.C:105
wordList toc() const
Return the table of contents.
Definition: dictionary.C:1018
IOoject and searching on triSurface.
Templated 2D Barycentric derived from VectorSpace. Has 3 components, one of which is redundant...
Definition: Barycentric2D.H:50
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
const pointField & points
A class for handling words, derived from string.
Definition: word.H:59
Container for searchableSurfaces.
Triangle with additional region number.
Definition: labelledTri.H:57
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
bool set(const Key &, const T &newElmt)
Assign a new hashedEntry, overwriting existing entries.
Definition: HashTableI.H:91
static const boundBox invertedBox
A very large inverted boundBox: min/max == +/- vGreat.
Definition: boundBox.H:82
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
const wordList & names() const
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:84
messageStream Info
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
A keyword and a list of tokens is an &#39;entry&#39;.
Definition: entry.H:65
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812
IOerror FatalIOError