layerParameters.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-2019 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 "layerParameters.H"
27 #include "polyBoundaryMesh.H"
28 #include "unitConversion.H"
29 #include "refinementSurfaces.H"
30 #include "searchableSurfaces.H"
31 #include "medialAxisMeshMover.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 Foam::scalar Foam::layerParameters::layerExpansionRatio
41 (
42  const label n,
43  const scalar totalOverFirst
44 ) const
45 {
46  if (n <= 1)
47  {
48  return 1;
49  }
50 
51  const scalar tol = 1e-8;
52 
53  if (mag(n - totalOverFirst) < tol)
54  {
55  return 1;
56  }
57 
58  const label maxIters = 100;
59 
60  // Calculate the bounds of the solution
61  scalar minR;
62  scalar maxR;
63 
64  if (totalOverFirst < n)
65  {
66  minR = 0;
67  maxR = pow(totalOverFirst/n, 1/(n-1));
68  }
69  else
70  {
71  minR = pow(totalOverFirst/n, 1/(n-1));
72  maxR = totalOverFirst/(n - 1);
73  }
74 
75  // Starting guess
76  scalar r = 0.5*(minR + maxR);
77 
78  for (label i = 0; i < maxIters; ++i)
79  {
80  const scalar prevr = r;
81  const scalar fx = pow(r, n) - totalOverFirst*r - (1 - totalOverFirst);
82  const scalar dfx = n*pow(r, n - 1) - totalOverFirst;
83 
84  r -= fx/dfx;
85 
86  if (mag(r - prevr) < tol)
87  {
88  break;
89  }
90  }
91 
92  return r;
93 }
94 
95 
96 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
97 
99 (
100  const dictionary& dict,
101  const polyBoundaryMesh& boundaryMesh
102 )
103 :
104  dict_(dict),
105  numLayers_(boundaryMesh.size(), -1),
106  relativeSizes_(dict.lookup("relativeSizes")),
107  layerSpec_(ILLEGAL),
108  firstLayerThickness_(boundaryMesh.size(), -123),
109  finalLayerThickness_(boundaryMesh.size(), -123),
110  thickness_(boundaryMesh.size(), -123),
111  expansionRatio_(boundaryMesh.size(), -123),
112  minThickness_
113  (
114  boundaryMesh.size(),
115  dict.lookup<scalar>("minThickness")
116  ),
117  mergeFaces_
118  (
119  boundaryMesh.size(),
120  dict.found("mergeFaces")
121  ? (
122  dict.lookup<bool>("mergeFaces")
123  ? mergeFace::yes
124  : mergeFace::no
125  )
126  : mergeFace::ifOnMeshedPatch
127  ),
128  featureAngle_(dict.lookup<scalar>("featureAngle")),
129  concaveAngle_
130  (
131  dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
132  ),
133  nGrow_(dict.lookup<label>("nGrow")),
134  maxFaceThicknessRatio_
135  (
136  dict.lookup<scalar>("maxFaceThicknessRatio")
137  ),
138  nBufferCellsNoExtrude_
139  (
140  dict.lookup<label>("nBufferCellsNoExtrude")
141  ),
142  nLayerIter_(dict.lookup<label>("nLayerIter")),
143  nRelaxedIter_(labelMax),
144  additionalReporting_(dict.lookupOrDefault("additionalReporting", false)),
145  meshShrinker_
146  (
147  dict.lookupOrDefault
148  (
149  "meshShrinker",
150  medialAxisMeshMover::typeName
151  )
152  )
153 {
154  // Detect layer specification mode
155 
156  label nSpec = 0;
157 
158  bool haveFirst = dict.found("firstLayerThickness");
159  if (haveFirst)
160  {
161  firstLayerThickness_ = scalarField
162  (
163  boundaryMesh.size(),
164  dict.lookup<scalar>("firstLayerThickness")
165  );
166  nSpec++;
167  }
168  bool haveFinal = dict.found("finalLayerThickness");
169  if (haveFinal)
170  {
171  finalLayerThickness_ = scalarField
172  (
173  boundaryMesh.size(),
174  dict.lookup<scalar>("finalLayerThickness")
175  );
176  nSpec++;
177  }
178  bool haveTotal = dict.found("thickness");
179  if (haveTotal)
180  {
181  thickness_ = scalarField
182  (
183  boundaryMesh.size(),
184  dict.lookup<scalar>("thickness")
185  );
186  nSpec++;
187  }
188  bool haveExp = dict.found("expansionRatio");
189  if (haveExp)
190  {
191  expansionRatio_ = scalarField
192  (
193  boundaryMesh.size(),
194  dict.lookup<scalar>("expansionRatio")
195  );
196  nSpec++;
197  }
198 
199 
200  if (haveFirst && haveTotal)
201  {
202  layerSpec_ = FIRST_AND_TOTAL;
203  Info<< "Layer thickness specified as first layer and overall thickness."
204  << endl;
205  }
206  else if (haveFirst && haveExp)
207  {
208  layerSpec_ = FIRST_AND_EXPANSION;
209  Info<< "Layer thickness specified as first layer and expansion ratio."
210  << endl;
211  }
212  else if (haveFinal && haveTotal)
213  {
214  layerSpec_ = FINAL_AND_TOTAL;
215  Info<< "Layer thickness specified as final layer and overall thickness."
216  << endl;
217  }
218  else if (haveFinal && haveExp)
219  {
220  layerSpec_ = FINAL_AND_EXPANSION;
221  Info<< "Layer thickness specified as final layer and expansion ratio."
222  << endl;
223  }
224  else if (haveTotal && haveExp)
225  {
226  layerSpec_ = TOTAL_AND_EXPANSION;
227  Info<< "Layer thickness specified as overall thickness"
228  << " and expansion ratio." << endl;
229  }
230 
231 
232  if (layerSpec_ == ILLEGAL || nSpec != 2)
233  {
235  (
236  dict
237  ) << "Over- or underspecified layer thickness."
238  << " Please specify" << nl
239  << " first layer thickness ('firstLayerThickness')"
240  << " and overall thickness ('thickness') or" << nl
241  << " first layer thickness ('firstLayerThickness')"
242  << " and expansion ratio ('expansionRatio') or" << nl
243  << " final layer thickness ('finalLayerThickness')"
244  << " and expansion ratio ('expansionRatio') or" << nl
245  << " final layer thickness ('finalLayerThickness')"
246  << " and overall thickness ('thickness') or" << nl
247  << " overall thickness ('thickness')"
248  << " and expansion ratio ('expansionRatio'"
249  << exit(FatalIOError);
250  }
251 
252 
253  dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
254 
255  if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
256  {
258  << "Layer iterations should be >= 0." << endl
259  << "nLayerIter:" << nLayerIter_
260  << " nRelaxedIter:" << nRelaxedIter_
261  << exit(FatalIOError);
262  }
263 
264 
265  const dictionary& layersDict = dict.subDict("layers");
266 
267  forAllConstIter(dictionary, layersDict, iter)
268  {
269  if (iter().isDict())
270  {
271  const keyType& key = iter().keyword();
272  const labelHashSet patchIDs
273  (
274  boundaryMesh.patchSet(List<wordRe>(1, wordRe(key)))
275  );
276 
277  if (patchIDs.size() == 0)
278  {
279  IOWarningInFunction(layersDict)
280  << "Layer specification for " << key
281  << " does not match any patch." << endl
282  << "Valid patches are " << boundaryMesh.names() << endl;
283  }
284  else
285  {
286  const dictionary& layerDict = iter().dict();
287 
288  forAllConstIter(labelHashSet, patchIDs, patchiter)
289  {
290  const label patchi = patchiter.key();
291 
292  numLayers_[patchi] =
293  layerDict.lookup<label>("nSurfaceLayers");
294 
295  switch (layerSpec_)
296  {
297  case FIRST_AND_TOTAL:
298  layerDict.readIfPresent
299  (
300  "firstLayerThickness",
301  firstLayerThickness_[patchi]
302  );
303  layerDict.readIfPresent
304  (
305  "thickness",
306  thickness_[patchi]
307  );
308  break;
309 
310  case FIRST_AND_EXPANSION:
311  layerDict.readIfPresent
312  (
313  "firstLayerThickness",
314  firstLayerThickness_[patchi]
315  );
316  layerDict.readIfPresent
317  (
318  "expansionRatio",
319  expansionRatio_[patchi]
320  );
321  break;
322 
323  case FINAL_AND_TOTAL:
324  layerDict.readIfPresent
325  (
326  "finalLayerThickness",
327  finalLayerThickness_[patchi]
328  );
329  layerDict.readIfPresent
330  (
331  "thickness",
332  thickness_[patchi]
333  );
334  break;
335 
336  case FINAL_AND_EXPANSION:
337  layerDict.readIfPresent
338  (
339  "finalLayerThickness",
340  finalLayerThickness_[patchi]
341  );
342  layerDict.readIfPresent
343  (
344  "expansionRatio",
345  expansionRatio_[patchi]
346  );
347  break;
348 
349  case TOTAL_AND_EXPANSION:
350  layerDict.readIfPresent
351  (
352  "thickness",
353  thickness_[patchi]
354  );
355  layerDict.readIfPresent
356  (
357  "expansionRatio",
358  expansionRatio_[patchi]
359  );
360  break;
361 
362  default:
364  (
365  dict
366  ) << "problem." << exit(FatalIOError);
367  break;
368  }
369 
370  layerDict.readIfPresent
371  (
372  "minThickness",
373  minThickness_[patchi]
374  );
375 
376  if (layerDict.found("mergeFaces"))
377  {
378  mergeFaces_[patchi] =
379  layerDict.lookup<bool>("mergeFaces")
381  : mergeFace::no;
382  }
383  }
384  }
385  }
386  }
387 }
388 
389 
390 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
391 
393 (
394  const label nLayers,
395  const scalar firstLayerThickness,
396  const scalar finalLayerThickness,
397  const scalar totalThickness,
398  const scalar expansionRatio
399 ) const
400 {
401  switch (layerSpec_)
402  {
403  case FIRST_AND_TOTAL:
404  case FINAL_AND_TOTAL:
405  case TOTAL_AND_EXPANSION:
406  {
407  return totalThickness;
408  }
409  break;
410 
411  case FIRST_AND_EXPANSION:
412  {
413  if (mag(expansionRatio-1) < small)
414  {
415  return firstLayerThickness * nLayers;
416  }
417  else
418  {
419  return firstLayerThickness
420  *(1 - pow(expansionRatio, nLayers))
421  /(1 - expansionRatio);
422  }
423  }
424  break;
425 
426  case FINAL_AND_EXPANSION:
427  {
428  if (mag(expansionRatio-1) < small)
429  {
430  return finalLayerThickness * nLayers;
431  }
432  else
433  {
434  const scalar invExpansion = 1.0/expansionRatio;
435 
436  return finalLayerThickness
437  *(1 - pow(invExpansion, nLayers))
438  /(1 - invExpansion);
439  }
440  }
441  break;
442 
443  default:
444  {
446  << exit(FatalError);
447  return -vGreat;
448  }
449  }
450 }
451 
452 
453 Foam::scalar Foam::layerParameters::layerExpansionRatio
454 (
455  const label nLayers,
456  const scalar firstLayerThickness,
457  const scalar finalLayerThickness,
458  const scalar totalThickness,
459  const scalar expansionRatio
460 ) const
461 {
462  switch (layerSpec_)
463  {
464  case FIRST_AND_EXPANSION:
465  case FINAL_AND_EXPANSION:
466  case TOTAL_AND_EXPANSION:
467  {
468  return expansionRatio;
469  }
470  break;
471 
472  case FIRST_AND_TOTAL:
473  {
474  return layerExpansionRatio
475  (
476  nLayers,
477  totalThickness/firstLayerThickness
478  );
479  }
480  break;
481 
482  case FINAL_AND_TOTAL:
483  {
484  return
485  1.0
486  /layerExpansionRatio
487  (
488  nLayers,
489  totalThickness/finalLayerThickness
490  );
491  }
492  break;
493 
494  default:
495  {
497  << "Illegal thickness specification" << exit(FatalError);
498  return -vGreat;
499  }
500  }
501 }
502 
503 
505 (
506  const label nLayers,
507  const scalar firstLayerThickness,
508  const scalar finalLayerThickness,
509  const scalar totalThickness,
510  const scalar expansionRatio
511 ) const
512 {
513  switch (layerSpec_)
514  {
515  case FIRST_AND_EXPANSION:
516  case FIRST_AND_TOTAL:
517  {
518  return firstLayerThickness;
519  }
520 
521  case FINAL_AND_EXPANSION:
522  {
523  return finalLayerThickness*pow(1.0/expansionRatio, nLayers-1);
524  }
525  break;
526 
527  case FINAL_AND_TOTAL:
528  {
529  const scalar r = layerExpansionRatio
530  (
531  nLayers,
532  firstLayerThickness,
533  finalLayerThickness,
534  totalThickness,
535  expansionRatio
536  );
537 
538  return finalLayerThickness/pow(r, nLayers-1);
539  }
540  break;
541 
542  case TOTAL_AND_EXPANSION:
543  {
544  const scalar r = finalLayerThicknessRatio
545  (
546  nLayers,
547  expansionRatio
548  );
549 
550  const scalar finalThickness = r*totalThickness;
551 
552  return finalThickness/pow(expansionRatio, nLayers-1);
553  }
554  break;
555 
556  default:
557  {
559  << "Illegal thickness specification" << exit(FatalError);
560  return -vGreat;
561  }
562  }
563 }
564 
565 
567 (
568  const label nLayers,
569  const scalar expansionRatio
570 ) const
571 {
572  if (nLayers > 0)
573  {
574  if (mag(expansionRatio-1) < small)
575  {
576  return 1.0/nLayers;
577  }
578  else
579  {
580  return
581  pow(expansionRatio, nLayers - 1)
582  *(1 - expansionRatio)
583  /(1 - pow(expansionRatio, nLayers));
584  }
585  }
586  else
587  {
588  return 0;
589  }
590 }
591 
592 
593 // ************************************************************************* //
bool found
label n
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:860
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:998
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:659
A class for handling keywords in dictionaries.
Definition: keyType.H:69
layerParameters(const dictionary &dict, const polyBoundaryMesh &)
Construct from dictionary.
scalar finalLayerThicknessRatio(const label nLayers, const scalar expansionRatio) const
Determine ratio of final layer thickness to.
const dictionary & dict() const
mergeFace
Enumeration defining whether to merge faces on a given patch. Read.
const scalarField & firstLayerThickness() const
Wanted thickness of the layer nearest to the wall.
scalar layerThickness(const label nLayers, const scalar firstLayerThickness, const scalar finalLayerThickness, const scalar totalThickness, const scalar expansionRatio) const
Determine overall thickness. Uses two of the four parameters.
Mesh motion solver that uses a medial axis algorithm to work out a fraction between the (nearest poin...
Foam::polyBoundaryMesh.
labelHashSet patchSet(const UList< wordRe > &patchNames, const bool warnNotFound=true, const bool usePatchGroups=true) const
Return the set of patch IDs corresponding to the given names.
wordList names() const
Return a list of patch names.
A wordRe is a word, but can also have a regular expression for matching words.
Definition: wordRe.H:77
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label patchi
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const doubleScalar e
Definition: doubleScalar.H:105
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
messageStream Info
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
IOerror FatalIOError
error FatalError
static const label labelMax
Definition: label.H:62
static const char nl
Definition: Ostream.H:260
dictionary dict
Unit conversion functions.