automatic.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) 2012-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 "automatic.H"
28 #include "triSurfaceMesh.H"
29 #include "vtkSurfaceWriter.H"
31 #include "Time.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(automatic, 0);
38  addToRunTimeSelectionTable(cellSizeCalculationType, automatic, dictionary);
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 void Foam::automatic::smoothField(triSurfaceScalarField& surf)
45 {
46  label nSmoothingIterations = 10;
47 
48  for (label iter = 0; iter < nSmoothingIterations; ++iter)
49  {
50  const pointField& faceCentres = surface_.faceCentres();
51 
52  forAll(surf, sI)
53  {
54  const labelList& faceFaces = surface_.faceFaces()[sI];
55 
56  const point& fC = faceCentres[sI];
57  const scalar value = surf[sI];
58 
59  scalar newValue = 0;
60  scalar totalDist = 0;
61 
62  label nFaces = 0;
63 
64  forAll(faceFaces, fI)
65  {
66  const label faceLabel = faceFaces[fI];
67  const point& faceCentre = faceCentres[faceLabel];
68 
69  const scalar faceValue = surf[faceLabel];
70  const scalar distance = mag(faceCentre - fC);
71 
72  newValue += faceValue/(distance + SMALL);
73 
74  totalDist += 1.0/(distance + SMALL);
75 
76  if (value < faceValue)
77  {
78  nFaces++;
79  }
80  }
81 
82  // Do not smooth out the peak values
83  if (nFaces == faceFaces.size())
84  {
85  continue;
86  }
87 
88  surf[sI] = newValue/totalDist;
89  }
90  }
91 }
92 
93 
94 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
95 
97 (
98  const dictionary& cellSizeCalcTypeDict,
99  const triSurfaceMesh& surface,
100  const scalar& defaultCellSize
101 )
102 :
103  cellSizeCalculationType
104  (
105  typeName,
106  cellSizeCalcTypeDict,
107  surface,
108  defaultCellSize
109  ),
110  coeffsDict_(cellSizeCalcTypeDict.subDict(typeName + "Coeffs")),
111  surfaceName_(surface.searchableSurface::name()),
112  readCurvature_(Switch(coeffsDict_.lookup("curvature"))),
113  curvatureFile_(coeffsDict_.lookup("curvatureFile")),
114  readFeatureProximity_(Switch(coeffsDict_.lookup("featureProximity"))),
115  featureProximityFile_(coeffsDict_.lookup("featureProximityFile")),
116  readInternalCloseness_(Switch(coeffsDict_.lookup("internalCloseness"))),
117  internalClosenessFile_(coeffsDict_.lookup("internalClosenessFile")),
118  curvatureCellSizeCoeff_
119  (
120  readScalar(coeffsDict_.lookup("curvatureCellSizeCoeff"))
121  ),
122  maximumCellSize_
123  (
124  readScalar(coeffsDict_.lookup("maximumCellSizeCoeff"))*defaultCellSize
125  )
126 {}
127 
128 
129 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
130 
132 {
133  Info<< indent
134  << "Calculating cell size on surface: " << surfaceName_ << endl;
135 
136  tmp<triSurfacePointScalarField> tPointCellSize
137  (
139  (
140  IOobject
141  (
142  surfaceName_ + ".cellSize",
143  surface_.searchableSurface::time().constant(),
144  "triSurface",
145  surface_.searchableSurface::time(),
148  ),
149  surface_,
150  dimLength,
151  scalarField(surface_.nPoints(), maximumCellSize_)
152  )
153  );
154 
155  triSurfacePointScalarField& pointCellSize = tPointCellSize.ref();
156 
157  if (readCurvature_)
158  {
159  Info<< indent
160  << "Reading curvature : " << curvatureFile_ << endl;
161 
163  (
164  IOobject
165  (
166  curvatureFile_,
167  surface_.searchableSurface::time().constant(),
168  "triSurface",
169  surface_.searchableSurface::time(),
172  ),
173  surface_,
174  dimLength,
175  true
176  );
177 
178  forAll(pointCellSize, pI)
179  {
180  pointCellSize[pI] =
181  min
182  (
183  1.0
184  /max
185  (
186  (1.0/curvatureCellSizeCoeff_)*mag(curvature[pI]),
187  1.0/maximumCellSize_
188  ),
189  pointCellSize[pI]
190  );
191  }
192  }
193 
194  PrimitivePatchInterpolation
195  <
196  PrimitivePatch<labelledTri, ::Foam::List, pointField, point>
197  > patchInterpolate(surface_);
198 
199  const Map<label>& meshPointMap = surface_.meshPointMap();
200 
201  if (readInternalCloseness_)
202  {
203  Info<< indent
204  << "Reading internal closeness: " << internalClosenessFile_ << endl;
205 
206  triSurfaceScalarField internalCloseness
207  (
208  IOobject
209  (
210  internalClosenessFile_,
211  surface_.searchableSurface::time().constant(),
212  "triSurface",
213  surface_.searchableSurface::time(),
216  ),
217  surface_,
218  dimLength,
219  true
220  );
221 
222  scalarField internalClosenessPointField
223  (
224  patchInterpolate.faceToPointInterpolate(internalCloseness)
225  );
226 
227  forAll(pointCellSize, pI)
228  {
229  pointCellSize[pI] =
230  min
231  (
232  internalClosenessPointField[meshPointMap[pI]],
233  pointCellSize[pI]
234  );
235  }
236  }
237 
238  if (readFeatureProximity_)
239  {
240  Info<< indent
241  << "Reading feature proximity : " << featureProximityFile_ << endl;
242 
243  triSurfaceScalarField featureProximity
244  (
245  IOobject
246  (
247  featureProximityFile_,
248  surface_.searchableSurface::time().constant(),
249  "triSurface",
250  surface_.searchableSurface::time(),
253  ),
254  surface_,
255  dimLength,
256  true
257  );
258 
259  scalarField featureProximityPointField
260  (
261  patchInterpolate.faceToPointInterpolate(featureProximity)
262  );
263 
264  forAll(pointCellSize, pI)
265  {
266  pointCellSize[pI] =
267  min
268  (
269  featureProximityPointField[meshPointMap[pI]],
270  pointCellSize[pI]
271  );
272  }
273  }
274 
275  //smoothField(surfaceCellSize);
276 
277  pointCellSize.write();
278 
279  if (debug)
280  {
281  faceList faces(surface_.size());
282 
283  forAll(surface_, fI)
284  {
285  faces[fI] = surface_.triSurface::operator[](fI).triFaceFace();
286  }
287 
288  vtkSurfaceWriter().write
289  (
290  surface_.searchableSurface::time().constant()/"triSurface",
291  surfaceName_.lessExt().name(),
292  surface_.points(),
293  faces,
294  "cellSize",
295  pointCellSize,
296  true,
297  true
298  );
299  }
300 
301  return tPointCellSize;
302 }
303 
304 
305 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Macros for easy insertion into run-time selection tables.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual tmp< triSurfacePointScalarField > load()
Load the cell size field.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
virtual tmp< pointField > points() const =0
Get the points that define the surface.
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
automatic(const dictionary &cellSizeCalcTypeDict, const triSurfaceMesh &surface, const scalar &defaultCellSize)
Construct from components.
vector point
Point is a vector.
Definition: point.H:41
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
virtual bool write() const
Write using setting from DB.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::DimensionedField< scalar, triSurfacePointGeoMesh > triSurfacePointScalarField
virtual label size() const =0
Range of local indices that can be returned.
A class for managing temporary objects.
Definition: PtrList.H:54
const searchableSurface & surface_
Reference to the searchableSurface that cellSizeFunction.
Foam::DimensionedField< scalar, triSurfaceGeoMesh > triSurfaceScalarField
const word & name() const
Return name.
Definition: IOobject.H:260
Namespace for OpenFOAM.