refinementLevel.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 refinementLevel
29
30Group
31 grpMeshAdvancedUtilities
32
33Description
34 Attempt to determine the refinement levels of a refined cartesian mesh.
35 Run BEFORE snapping.
36
37 Writes
38 - volScalarField 'refinementLevel' with current refinement level.
39 - cellSet 'refCells' which are the cells that need to be refined to satisfy
40 2:1 refinement.
41
42 Works by dividing cells into volume bins.
43
44\*---------------------------------------------------------------------------*/
45
46#include "argList.H"
47#include "Time.H"
48#include "polyMesh.H"
49#include "cellSet.H"
50#include "SortableList.H"
51#include "labelIOList.H"
52#include "fvMesh.H"
53#include "volFields.H"
54
55using namespace Foam;
56
57// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58
59// Return true if any cells had to be split to keep a difference between
60// neighbouring refinement levels < limitDiff. Puts cells into refCells and
61// update refLevel to account for refinement.
62bool limitRefinementLevel
63(
64 const primitiveMesh& mesh,
65 labelList& refLevel,
66 cellSet& refCells
67)
68{
69 const labelListList& cellCells = mesh.cellCells();
70
71 label oldNCells = refCells.size();
72
73 forAll(cellCells, celli)
74 {
75 const labelList& cCells = cellCells[celli];
76
77 forAll(cCells, i)
78 {
79 if (refLevel[cCells[i]] > (refLevel[celli]+1))
80 {
81 // Found neighbour with >=2 difference in refLevel.
82 refCells.insert(celli);
83 refLevel[celli]++;
84 break;
85 }
86 }
87 }
88
89 if (refCells.size() > oldNCells)
90 {
91 Info<< "Added an additional " << refCells.size() - oldNCells
92 << " cells to satisfy 1:2 refinement level"
93 << endl;
94
95 return true;
96 }
97
98 return false;
99}
100
101
102int main(int argc, char *argv[])
103{
104 argList::addNote
105 (
106 "Attempt to determine refinement levels of a refined cartesian mesh.\n"
107 "Run BEFORE snapping!"
108 );
109
110 argList::addBoolOption
111 (
112 "readLevel",
113 "Read level from refinementLevel file"
114 );
115
116 #include "setRootCase.H"
117 #include "createTime.H"
118 #include "createPolyMesh.H"
119
120 Info<< "Dividing cells into bins depending on cell volume.\nThis will"
121 << " correspond to refinement levels for a mesh with only 2x2x2"
122 << " refinement\n"
123 << "The upper range for every bin is always 1.1 times the lower range"
124 << " to allow for some truncation error."
125 << nl << endl;
126
127 const bool readLevel = args.found("readLevel");
128
129 const scalarField& vols = mesh.cellVolumes();
130
131 SortableList<scalar> sortedVols(vols);
132
133 // All cell labels, sorted per bin.
135
136 // Lower/upper limits
137 DynamicList<scalar> lowerLimits;
138 DynamicList<scalar> upperLimits;
139
140 // Create bin0. Have upperlimit as factor times lowerlimit.
142 lowerLimits.append(sortedVols[0]);
143 upperLimits.append(1.1 * lowerLimits.last());
144
145 forAll(sortedVols, i)
146 {
147 if (sortedVols[i] > upperLimits.last())
148 {
149 // New value outside of current bin
150
151 // Shrink old bin.
152 DynamicList<label>& bin = bins.last();
153
154 bin.shrink();
155
156 Info<< "Collected " << bin.size() << " elements in bin "
157 << lowerLimits.last() << " .. "
158 << upperLimits.last() << endl;
159
160 // Create new bin.
162 lowerLimits.append(sortedVols[i]);
163 upperLimits.append(1.1 * lowerLimits.last());
164
165 Info<< "Creating new bin " << lowerLimits.last()
166 << " .. " << upperLimits.last()
167 << endl;
168 }
169
170 // Append to current bin.
171 DynamicList<label>& bin = bins.last();
172
173 bin.append(sortedVols.indices()[i]);
174 }
175 Info<< endl;
176
177 bins.last().shrink();
178 bins.shrink();
179 lowerLimits.shrink();
180 upperLimits.shrink();
181
182
183 //
184 // Write to cellSets.
185 //
186
187 Info<< "Volume bins:" << nl;
188 forAll(bins, binI)
189 {
190 const DynamicList<label>& bin = bins[binI];
191
192 cellSet cells(mesh, "vol" + name(binI), bin.size());
193 cells.insert(bin);
194
195 Info<< " " << lowerLimits[binI] << " .. " << upperLimits[binI]
196 << " : writing " << bin.size() << " cells to cellSet "
197 << cells.name() << endl;
198
199 cells.write();
200 }
201
202
203
204 //
205 // Convert bins into refinement level.
206 //
207
208
209 // Construct fvMesh to be able to construct volScalarField
210
211 fvMesh fMesh
212 (
214 (
215 polyMesh::defaultRegion,
216 runTime.timeName(),
217 runTime,
218 IOobject::NO_READ,
219 IOobject::AUTO_WRITE
220 ),
221 pointField(mesh.points()), // Could we safely re-use the data?
222 faceList(mesh.faces()),
223 cellList(mesh.cells())
224 );
225
226 // Add the boundary patches
227 const polyBoundaryMesh& patches = mesh.boundaryMesh();
228
229 List<polyPatch*> p(patches.size());
230
231 forAll(p, patchi)
232 {
233 p[patchi] = patches[patchi].clone(fMesh.boundaryMesh()).ptr();
234 }
235
236 fMesh.addFvPatches(p);
237
238
239 // Refinement level
240 IOobject refHeader
241 (
242 "refinementLevel",
243 runTime.timeName(),
244 polyMesh::defaultRegion,
245 runTime
246 );
247
248 if (!readLevel && refHeader.typeHeaderOk<labelIOList>(true))
249 {
251 << "Detected " << refHeader.name() << " file in "
252 << polyMesh::defaultRegion << " directory. Please remove to"
253 << " recreate it or use the -readLevel option to use it"
254 << endl;
255 return 1;
256 }
257
258
259 labelIOList refLevel
260 (
262 (
263 "refinementLevel",
264 runTime.timeName(),
265 mesh,
266 IOobject::NO_READ,
267 IOobject::AUTO_WRITE
268 ),
269 labelList(mesh.nCells(), Zero)
270 );
271
272 if (readLevel)
273 {
274 refLevel = labelIOList(refHeader);
275 }
276
277 // Construct volScalarField with same info for post processing
278 volScalarField postRefLevel
279 (
281 (
282 "refinementLevel",
283 runTime.timeName(),
284 mesh,
285 IOobject::NO_READ,
286 IOobject::NO_WRITE
287 ),
288 fMesh,
289 dimensionedScalar(dimless/dimTime, Zero)
290 );
291
292 // Set cell values
293 forAll(bins, binI)
294 {
295 const DynamicList<label>& bin = bins[binI];
296
297 forAll(bin, i)
298 {
299 refLevel[bin[i]] = bins.size() - binI - 1;
300 postRefLevel[bin[i]] = refLevel[bin[i]];
301 }
302 }
303
304 volScalarField::Boundary& postRefLevelBf =
305 postRefLevel.boundaryFieldRef();
306
307 // For volScalarField: set boundary values to same as cell.
308 // Note: could also put
309 // zeroGradient b.c. on postRefLevel and do evaluate.
310 forAll(postRefLevel.boundaryField(), patchi)
311 {
312 const polyPatch& pp = patches[patchi];
313
314 fvPatchScalarField& bField = postRefLevelBf[patchi];
315
316 Info<< "Setting field for patch "<< endl;
317
318 forAll(bField, facei)
319 {
320 label own = mesh.faceOwner()[pp.start() + facei];
321
322 bField[facei] = postRefLevel[own];
323 }
324 }
325
326 Info<< "Determined current refinement level and writing to "
327 << postRefLevel.name() << " (as volScalarField; for post processing)"
328 << nl
329 << polyMesh::defaultRegion/refLevel.name()
330 << " (as labelIOList; for meshing)" << nl
331 << endl;
332
333 refLevel.write();
334 postRefLevel.write();
335
336
337 // Find out cells to refine to keep to 2:1 refinement level restriction
338
339 // Cells to refine
340 cellSet refCells(mesh, "refCells", 100);
341
342 while
343 (
344 limitRefinementLevel
345 (
346 mesh,
347 refLevel, // current refinement level
348 refCells // cells to refine
349 )
350 )
351 {}
352
353 if (refCells.size())
354 {
355 Info<< "Collected " << refCells.size() << " cells that need to be"
356 << " refined to get closer to overall 2:1 refinement level limit"
357 << nl
358 << "Written cells to be refined to cellSet " << refCells.name()
359 << nl << endl;
360
361 refCells.write();
362
363 Info<< "After refinement this tool can be run again to see if the 2:1"
364 << " limit is observed all over the mesh" << nl << endl;
365 }
366 else
367 {
368 Info<< "All cells in the mesh observe the 2:1 refinement level limit"
369 << nl << endl;
370 }
371
372 Info<< "\nEnd\n" << endl;
373 return 0;
374}
375
376
377// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:63
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
A collection of cell labels.
Definition: cellSet.H:54
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
virtual bool write(const bool valid=true) const
Write using setting from DB.
volScalarField & p
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
const cellShapeList & cells
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333