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