refinementRegions.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-2021 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 "refinementRegions.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<refinementRegions::refineMode, 5>::
40 names[] =
41 {
42  "inside",
43  "outside",
44  "distance",
45  "insideSpan",
46  "outsideSpan"
47 };
48 
49 }
50 
52  Foam::refinementRegions::refineModeNames_;
53 
54 
55 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 
57 void Foam::refinementRegions::setAndCheckLevels
58 (
59  const label shelli,
60  const dictionary& dict
61 )
62 {
63  const searchableSurface& shell = allGeometry_[shells_[shelli]];
64 
65  if
66  (
67  modes_[shelli] == refineMode::inside
68  || modes_[shelli] == refineMode::outside
69  )
70  {
71  if (!allGeometry_[shells_[shelli]].hasVolumeType())
72  {
74  << "Shell " << shell.name()
75  << " is not closed so testing for '"
76  << refineModeNames_[modes_[shelli]]
77  << "' may fail." << endl;
78  }
79 
80  distances_[shelli].setSize(0);
81  levels_[shelli].setSize(1);
82 
83  if (dict.found("levels") && !(dict.found("level")))
84  {
85  // Support 'levels' for backward compatibility
86  const List<Tuple2<scalar, label>> distLevels(dict.lookup("levels"));
87 
88  if (distLevels.size() != 1)
89  {
91  << "For refinement mode "
92  << refineModeNames_[modes_[shelli]]
93  << " specify only one distance+level."
94  << " (its distance gets discarded)"
95  << exit(FatalError);
96  }
97 
98  levels_[shelli][0] = distLevels[0].second();
99  }
100  else
101  {
102  if (dict.found("levels"))
103  {
104  IOWarningInFunction(dict)
105  << "Found both 'level' and 'levels' entries, using 'level'."
106  << endl;
107  }
108 
109  levels_[shelli][0] = readLabel(dict.lookup("level"));
110  }
111 
112  if (modes_[shelli] == refineMode::inside)
113  {
114  Info<< "Refinement level " << levels_[shelli][0]
115  << " for all cells inside " << shell.name() << endl;
116  }
117  else
118  {
119  Info<< "Refinement level " << levels_[shelli][0]
120  << " for all cells outside " << shell.name() << endl;
121  }
122  }
123  else if
124  (
125  modes_[shelli] == refineMode::insideSpan
126  || modes_[shelli] == refineMode::outsideSpan
127  )
128  {
129  if (!allGeometry_[shells_[shelli]].hasVolumeType())
130  {
132  << "Shell " << shell.name()
133  << " is not closed so testing for '"
134  << refineModeNames_[modes_[shelli]]
135  << "' may fail." << endl;
136  }
137 
138  distances_[shelli].setSize(1);
139  levels_[shelli].setSize(1);
140  const Tuple2<scalar, label> distLevel(dict.lookup("level"));
141  distances_[shelli][0] = distLevel.first();
142  levels_[shelli][0] = distLevel.second();
143 
144  if (modes_[shelli] == refineMode::insideSpan)
145  {
146  Info<< "Refinement level " << levels_[shelli][0]
147  << " for all cells inside " << shell.name()
148  << " within distance " << distances_[shelli][0] << endl;
149  }
150  else
151  {
152  Info<< "Refinement level " << levels_[shelli][0]
153  << " for all cells outside " << shell.name()
154  << " within distance " << distances_[shelli][0] << endl;
155  }
156  }
157  else
158  {
159  const List<Tuple2<scalar, label>> distLevels(dict.lookup("levels"));
160 
161  // Extract information into separate distance and level
162  distances_[shelli].setSize(distLevels.size());
163  levels_[shelli].setSize(distLevels.size());
164 
165  forAll(distLevels, j)
166  {
167  distances_[shelli][j] = distLevels[j].first();
168  levels_[shelli][j] = distLevels[j].second();
169 
170  // Check in incremental order
171  if (j > 0)
172  {
173  if
174  (
175  (distances_[shelli][j] <= distances_[shelli][j-1])
176  || (levels_[shelli][j] > levels_[shelli][j-1])
177  )
178  {
180  << "For refinement mode "
181  << refineModeNames_[modes_[shelli]]
182  << " : Refinement should be specified in order"
183  << " of increasing distance"
184  << " (and decreasing refinement level)." << endl
185  << "Distance:" << distances_[shelli][j]
186  << " refinementLevel:" << levels_[shelli][j]
187  << exit(FatalError);
188  }
189  }
190  }
191 
192  if (modes_[shelli] == refineMode::distance)
193  {
194  Info<< "Refinement level according to distance to "
195  << shell.name() << endl;
196 
197  forAll(levels_[shelli], j)
198  {
199  Info<< " level " << levels_[shelli][j]
200  << " for all cells within " << distances_[shelli][j]
201  << " metre." << endl;
202  }
203  }
204  }
205 }
206 
207 
208 void Foam::refinementRegions::orient()
209 {
210  // Determine outside point.
211  boundBox overallBb = boundBox::invertedBox;
212 
213  bool hasSurface = false;
214 
215  forAll(shells_, shelli)
216  {
217  const searchableSurface& s = allGeometry_[shells_[shelli]];
218 
219  if (modes_[shelli] != refineMode::distance && isA<triSurfaceMesh>(s))
220  {
221  const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
222 
223  if (shell.triSurface::size())
224  {
225  const pointField& points = shell.points();
226 
227  hasSurface = true;
228 
229  boundBox shellBb(points[0], points[0]);
230  // Assume surface is compact!
231  forAll(points, i)
232  {
233  const point& pt = points[i];
234  shellBb.min() = min(shellBb.min(), pt);
235  shellBb.max() = max(shellBb.max(), pt);
236  }
237 
238  overallBb.min() = min(overallBb.min(), shellBb.min());
239  overallBb.max() = max(overallBb.max(), shellBb.max());
240  }
241  }
242  }
243 
244  if (hasSurface)
245  {
246  const point outsidePt = overallBb.max() + overallBb.span();
247 
248  // Info<< "Using point " << outsidePt << " to orient shells" << endl;
249 
250  forAll(shells_, shelli)
251  {
252  const searchableSurface& s = allGeometry_[shells_[shelli]];
253 
254  if
255  (
256  modes_[shelli] != refineMode::distance
257  && isA<triSurfaceMesh>(s)
258  )
259  {
260  triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
261  (
262  refCast<const triSurfaceMesh>(s)
263  );
264 
265  // Flip surface so outsidePt is outside.
266  bool anyFlipped = orientedSurface::orient
267  (
268  shell,
269  outsidePt,
270  true
271  );
272 
273  if (anyFlipped)
274  {
275  // orientedSurface will have done a clearOut of the surface.
276  // we could do a clearout of the triSurfaceMeshes::trees()
277  // but these aren't affected by orientation
278  // (except for cached
279  // sideness which should not be set at this point.
280  // !!Should check!)
281 
282  Info<< "refinementRegions : Flipped orientation of surface "
283  << s.name()
284  << " so point " << outsidePt << " is outside." << endl;
285  }
286  }
287  }
288  }
289 }
290 
291 
292 Foam::scalar Foam::refinementRegions::interpolate
293 (
294  const triSurfaceMesh& tsm,
295  const triSurfacePointScalarField& closeness,
296  const point& pt,
297  const label index
298 ) const
299 {
300  const barycentric2D bary
301  (
303  (
304  tsm.points(),
305  tsm.triSurface::operator[](index)
306  ).pointToBarycentric(pt)
307  );
308 
309  const labelledTri& lf = tsm.localFaces()[index];
310 
311  return closeness[lf[0]]*bary[0]
312  + closeness[lf[1]]*bary[1]
313  + closeness[lf[2]]*bary[2];
314 }
315 
316 
317 void Foam::refinementRegions::findHigherLevel
318 (
319  const pointField& pt,
320  const label shelli,
321  const scalar level0EdgeLength,
322  labelList& maxLevel
323 ) const
324 {
325  const labelList& levels = levels_[shelli];
326 
327  if (modes_[shelli] == refineMode::distance)
328  {
329  // Distance mode.
330 
331  const scalarField& distances = distances_[shelli];
332 
333  // Collect all those points that have a current maxLevel less than
334  // (any of) the shell. Also collect the furthest distance allowable
335  // to any shell with a higher level.
336 
337  pointField candidates(pt.size());
338  labelList candidateMap(pt.size());
339  scalarField candidateDistSqr(pt.size());
340  label candidatei = 0;
341 
342  forAll(maxLevel, pointi)
343  {
344  forAllReverse(levels, leveli)
345  {
346  if (levels[leveli] > maxLevel[pointi])
347  {
348  candidates[candidatei] = pt[pointi];
349  candidateMap[candidatei] = pointi;
350  candidateDistSqr[candidatei] = sqr(distances[leveli]);
351  candidatei++;
352  break;
353  }
354  }
355  }
356  candidates.setSize(candidatei);
357  candidateMap.setSize(candidatei);
358  candidateDistSqr.setSize(candidatei);
359 
360  // Do the expensive nearest test only for the candidate points.
362  allGeometry_[shells_[shelli]].findNearest
363  (
364  candidates,
365  candidateDistSqr,
366  nearInfo
367  );
368 
369  // Update maxLevel
370  forAll(nearInfo, candidatei)
371  {
372  if (nearInfo[candidatei].hit())
373  {
374  // Check which level it actually is in.
375  label minDistI = findLower
376  (
377  distances,
378  mag(nearInfo[candidatei].hitPoint()-candidates[candidatei])
379  );
380 
381  label pointi = candidateMap[candidatei];
382 
383  // pt is in between shell[minDistI] and shell[minDistI+1]
384  maxLevel[pointi] = levels[minDistI+1];
385  }
386  }
387  }
388  else if
389  (
390  modes_[shelli] == refineMode::insideSpan
391  || modes_[shelli] == refineMode::outsideSpan
392  )
393  {
394  const triSurfaceMesh& tsm =
395  refCast<const triSurfaceMesh>(allGeometry_[shells_[shelli]]);
396 
397  // Collect all those points that have a current maxLevel less than
398  // the maximum and furthest distance allowable for the shell.
399 
400  pointField candidates(pt.size());
401  labelList candidateMap(pt.size());
402  scalarField candidateDistSqr(pt.size());
403  label candidatei = 0;
404 
405  forAll(pt, pointi)
406  {
407  if (levels[0] > maxLevel[pointi])
408  {
409  candidates[candidatei] = pt[pointi];
410  candidateMap[candidatei] = pointi;
411  candidateDistSqr[candidatei] = sqr(distances_[shelli][0]);
412  candidatei++;
413  }
414  }
415  candidates.setSize(candidatei);
416  candidateMap.setSize(candidatei);
417  candidateDistSqr.setSize(candidatei);
418 
419  // Do the expensive nearest test only for the candidate points.
421  tsm.findNearest
422  (
423  candidates,
424  candidateDistSqr,
425  nearInfo
426  );
427 
428  // Minimum span for the maximum level specified
429  const scalar minSpan
430  (
431  cellsAcrossSpan_[shelli]*level0EdgeLength/(1 << (levels[0] - 1))
432  );
433 
434  // Update maxLevel
435  forAll(nearInfo, candidatei)
436  {
437  if (nearInfo[candidatei].hit())
438  {
439  const scalar span
440  (
442  (
443  tsm,
444  closeness_[shelli],
445  nearInfo[candidatei].rawPoint(),
446  nearInfo[candidatei].index()
447  )
448  );
449 
450  if (span > minSpan)
451  {
452  const label level
453  (
454  log2(cellsAcrossSpan_[shelli]*level0EdgeLength/span) + 1
455  );
456 
457  maxLevel[candidateMap[candidatei]] = min(levels[0], level);
458  }
459  else
460  {
461  maxLevel[candidateMap[candidatei]] = levels[0];
462  }
463  }
464  }
465  }
466  else
467  {
468  // Inside/outside mode
469 
470  // Collect all those points that have a current maxLevel less than the
471  // shell.
472 
473  pointField candidates(pt.size());
474  labelList candidateMap(pt.size());
475  label candidatei = 0;
476 
477  forAll(maxLevel, pointi)
478  {
479  if (levels[0] > maxLevel[pointi])
480  {
481  candidates[candidatei] = pt[pointi];
482  candidateMap[candidatei] = pointi;
483  candidatei++;
484  }
485  }
486  candidates.setSize(candidatei);
487  candidateMap.setSize(candidatei);
488 
489  // Do the expensive nearest test only for the candidate points.
490  List<volumeType> volType;
491  allGeometry_[shells_[shelli]].getVolumeType(candidates, volType);
492 
493  forAll(volType, i)
494  {
495  label pointi = candidateMap[i];
496 
497  if
498  (
499  (
500  modes_[shelli] == refineMode::inside
501  && volType[i] == volumeType::inside
502  )
503  || (
504  modes_[shelli] == refineMode::outside
505  && volType[i] == volumeType::outside
506  )
507  )
508  {
509  maxLevel[pointi] = levels[0];
510  }
511  }
512  }
513 }
514 
515 
516 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
517 
519 (
520  const searchableSurfaces& allGeometry,
521  const dictionary& shellsDict
522 )
523 :
524  allGeometry_(allGeometry)
525 {
526  // Wildcard specification : loop over all surfaces and try to find a match.
527 
528  // Count number of shells.
529  label shelli = 0;
530  forAll(allGeometry.names(), geomi)
531  {
532  const word& geomName = allGeometry_.names()[geomi];
533 
534  if (shellsDict.found(geomName))
535  {
536  shelli++;
537  }
538  }
539 
540  // Size lists
541  shells_.setSize(shelli);
542  modes_.setSize(shelli);
543  distances_.setSize(shelli);
544  levels_.setSize(shelli);
545  cellsAcrossSpan_.setSize(shelli);
546  closeness_.setSize(shelli);
547 
548  HashSet<word> unmatchedKeys(shellsDict.toc());
549  shelli = 0;
550 
551  forAll(allGeometry_.names(), geomi)
552  {
553  const word& geomName = allGeometry_.names()[geomi];
554 
555  const entry* ePtr = shellsDict.lookupEntryPtr(geomName, false, true);
556 
557  if (ePtr)
558  {
559  const dictionary& dict = ePtr->dict();
560  unmatchedKeys.erase(ePtr->keyword());
561 
562  shells_[shelli] = geomi;
563  modes_[shelli] = refineModeNames_.read(dict.lookup("mode"));
564 
565  // Read pairs of distance+level
566  setAndCheckLevels(shelli, dict);
567 
568  if
569  (
570  modes_[shelli] == refineMode::insideSpan
571  || modes_[shelli] == refineMode::outsideSpan
572  )
573  {
574  const searchableSurface& surface = allGeometry_[geomi];
575 
576  if (isA<triSurfaceMesh>(surface))
577  {
578  dict.lookup("cellsAcrossSpan") >> cellsAcrossSpan_[shelli];
579 
580  const triSurfaceMesh& tsm =
581  refCast<const triSurfaceMesh>(surface);
582 
583  closeness_.set
584  (
585  shelli,
587  (
588  IOobject
589  (
590  surface.name()
591  + ".closeness."
592  + (
593  modes_[shelli] == refineMode::insideSpan
594  ? "internal"
595  : "external"
596  )
597  + "PointCloseness",
598  surface.searchableSurface::time().constant(),
600  (
601  surface.searchableSurface::time()
602  ),
603  surface.searchableSurface::time(),
605  ),
606  tsm,
607  dimLength,
608  true
609  )
610  );
611  }
612  else
613  {
614  FatalIOErrorInFunction(shellsDict)
615  << "Surface " << surface.name()
616  << " is not a triSurface as required by"
617  " refinement modes "
618  << refineModeNames_[refineMode::insideSpan]
619  << " and " << refineModeNames_[refineMode::outsideSpan]
620  << exit(FatalIOError);
621  }
622  }
623 
624  shelli++;
625  }
626  }
627 
628  if (unmatchedKeys.size() > 0)
629  {
631  (
632  shellsDict
633  ) << "Not all entries in refinementRegions dictionary were used."
634  << " The following entries were not used : "
635  << unmatchedKeys.sortedToc()
636  << endl;
637  }
638 
639  // Orient shell surfaces before any searching is done. Note that this
640  // only needs to be done for inside or outside. Orienting surfaces
641  // constructs lots of addressing which we want to avoid.
642  orient();
643 }
644 
645 
646 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
647 
649 {
650  label overallMax = 0;
651 
652  forAll(levels_, shelli)
653  {
654  overallMax = max(overallMax, max(levels_[shelli]));
655  }
656 
657  return overallMax;
658 }
659 
660 
661 void Foam::refinementRegions::findHigherLevel
662 (
663  const pointField& pt,
664  const labelList& ptLevel,
665  const scalar level0EdgeLength,
666  labelList& maxLevel
667 ) const
668 {
669  // Maximum level of any shell. Start off with level of point.
670  maxLevel = ptLevel;
671 
672  forAll(shells_, shelli)
673  {
674  findHigherLevel(pt, shelli, level0EdgeLength, maxLevel);
675  }
676 }
677 
678 
679 // ************************************************************************* //
A HashTable with keys but without contents.
Definition: HashSet.H:59
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:663
#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
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const word & name() const
Return name.
Definition: IOobject.H:315
const entry * lookupEntryPtr(const word &, bool recursive, bool patternMatch) const
Find and return an entry data stream pointer if present.
Definition: dictionary.C:702
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:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
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
const Type1 & first() const
Return first.
Definition: Tuple2.H:116
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.
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
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
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:1115
const dimensionSet dimLength
A surface geometry formed of discrete facets, e.g. triangles and/or quadrilaterals, defined in a file using formats such as Wavefront OBJ, or stereolithography STL.
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))
refinementRegions(const searchableSurfaces &allGeometry, const dictionary &shellsDict)
Construct from geometry and dictionary.
const pointField & points
static const word & geometryDir()
Return the geometry directory name.
A class for handling words, derived from string.
Definition: word.H:59
Container for searchableSurfaces.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Triangle with additional region number.
Definition: labelledTri.H:57
label readLabel(Istream &is)
Definition: label.H:64
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum point defining 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
#define WarningInFunction
Report a warning using Foam::Warning.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
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 point defining the bounding box.
Definition: boundBoxI.H:54
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
label maxLevel() const
Highest shell level.
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:864
IOerror FatalIOError