multiDirRefinement.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) 2015-2020 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
27\*---------------------------------------------------------------------------*/
28
29#include "multiDirRefinement.H"
30#include "polyMesh.H"
31#include "Time.H"
32#include "undoableMeshCutter.H"
33#include "hexCellLooper.H"
34#include "geomCellLooper.H"
35#include "directions.H"
36#include "hexRef8.H"
37#include "mapPolyMesh.H"
38#include "polyTopoChange.H"
39#include "cellModel.H"
40
41// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42
43namespace Foam
44{
46}
47
48
49// * * * * * * * * * * * * * Private Statc Functions * * * * * * * * * * * * //
50
51// Update refCells pattern for split cells. Note refCells is current
52// list of cells to refine (these should all have been refined)
53void Foam::multiDirRefinement::addCells
54(
55 const Map<label>& splitMap,
56 List<refineCell>& refCells
57)
58{
59 label newRefI = refCells.size();
60
61 label oldSize = refCells.size();
62
63 refCells.setSize(newRefI + splitMap.size());
64
65 for (label refI = 0; refI < oldSize; refI++)
66 {
67 const refineCell& refCell = refCells[refI];
68
69 const auto iter = splitMap.cfind(refCell.cellNo());
70
71 if (!iter.found())
72 {
74 << "Problem : cannot find added cell for cell "
75 << refCell.cellNo() << endl
76 << abort(FatalError);
77 }
78
79 refCells[newRefI++] = refineCell(iter.val(), refCell.direction());
80 }
81}
82
83
84// Update vectorField for all the new cells. Takes over value of
85// original cell.
87(
88 const Map<label>& splitMap,
90)
91{
92 field.setSize(field.size() + splitMap.size());
93
94 forAllConstIters(splitMap, iter)
95 {
96 field[iter.val()] = field[iter.key()];
97 }
98}
99
100
101// Add added cells to labelList
102void Foam::multiDirRefinement::addCells
103(
104 const Map<label>& splitMap,
105 labelList& labels
106)
107{
108 label newCelli = labels.size();
109
110 labels.setSize(labels.size() + splitMap.size());
111
112 forAllConstIters(splitMap, iter)
113 {
114 labels[newCelli++] = iter.val();
115 }
116}
117
118
119// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
120
121void Foam::multiDirRefinement::addCells
122(
123 const primitiveMesh& mesh,
124 const Map<label>& splitMap
125)
126{
127 // Construct inverse addressing: from new to original cell.
128 labelList origCell(mesh.nCells(), -1);
129
130 forAll(addedCells_, celli)
131 {
132 const labelList& added = addedCells_[celli];
133
134 forAll(added, i)
135 {
136 label slave = added[i];
137
138 if (origCell[slave] == -1)
139 {
140 origCell[slave] = celli;
141 }
142 else if (origCell[slave] != celli)
143 {
145 << "Added cell " << slave << " has two different masters:"
146 << origCell[slave] << " , " << celli
147 << abort(FatalError);
148 }
149 }
150 }
151
152
153 forAllConstIters(splitMap, iter)
154 {
155 label masterI = iter.key();
156 const label newCelli = iter.val();
157
158 while (origCell[masterI] != -1 && origCell[masterI] != masterI)
159 {
160 masterI = origCell[masterI];
161 }
162
163 if (masterI >= addedCells_.size())
164 {
166 << "Map of added cells contains master cell " << masterI
167 << " which is not a valid cell number" << endl
168 << "This means that the mesh is not consistent with the"
169 << " done refinement" << endl
170 << "newCell:" << newCelli << abort(FatalError);
171 }
172
173 labelList& added = addedCells_[masterI];
174
175 if (added.empty())
176 {
177 added.setSize(2);
178 added[0] = masterI;
179 added[1] = newCelli;
180 }
181 else if (!added.found(newCelli))
182 {
183 const label sz = added.size();
184 added.setSize(sz + 1);
185 added[sz] = newCelli;
186 }
187 }
188}
189
190
191Foam::labelList Foam::multiDirRefinement::splitOffHex(const primitiveMesh& mesh)
192{
193 const cellModel& hex = cellModel::ref(cellModel::HEX);
194
196
197 // Split cellLabels_ into two lists.
198
199 labelList nonHexLabels(cellLabels_.size());
200 label nonHexI = 0;
201
202 labelList hexLabels(cellLabels_.size());
203 label hexI = 0;
204
205 forAll(cellLabels_, i)
206 {
207 label celli = cellLabels_[i];
208
209 if (cellShapes[celli].model() == hex)
210 {
211 hexLabels[hexI++] = celli;
212 }
213 else
214 {
215 nonHexLabels[nonHexI++] = celli;
216 }
217 }
218
219 nonHexLabels.setSize(nonHexI);
220
221 cellLabels_.transfer(nonHexLabels);
222
223 hexLabels.setSize(hexI);
224
225 return hexLabels;
226}
227
228
229void Foam::multiDirRefinement::refineHex8
230(
231 polyMesh& mesh,
232 const labelList& hexCells,
233 const bool writeMesh
234)
235{
236 if (debug)
237 {
238 Pout<< "multiDirRefinement : Refining hexes " << hexCells.size()
239 << endl;
240 }
241
242 hexRef8 hexRefiner
243 (
244 mesh,
245 labelList(mesh.nCells(), Zero), // cellLevel
246 labelList(mesh.nPoints(), Zero), // pointLevel
247 refinementHistory
248 (
249 IOobject
250 (
251 "refinementHistory",
254 mesh,
257 false
258 ),
259 List<refinementHistory::splitCell8>(0),
260 labelList(0),
261 false
262 ) // refinement history
263 );
264
265 polyTopoChange meshMod(mesh);
266
267 labelList consistentCells
268 (
269 hexRefiner.consistentRefinement
270 (
271 hexCells,
272 true // buffer layer
273 )
274 );
275
276 // Check that consistentRefinement hasn't added cells
277 {
278 // Create count 1 for original cells
279 Map<label> hexCellSet(2*hexCells.size());
280 forAll(hexCells, i)
281 {
282 hexCellSet.insert(hexCells[i], 1);
283 }
284
285 // Increment count
286 forAll(consistentCells, i)
287 {
288 const label celli = consistentCells[i];
289
290 auto iter = hexCellSet.find(celli);
291
292 if (iter.found())
293 {
294 iter.val() = 2;
295 }
296 else
297 {
299 << "Resulting mesh would not satisfy 2:1 ratio"
300 << " when refining cell " << celli << abort(FatalError);
301 }
302 }
303
304 // Check if all been visited (should always be since
305 // consistentRefinement set up to extend set.
306 forAllConstIters(hexCellSet, iter)
307 {
308 if (iter.val() != 2)
309 {
311 << "Resulting mesh would not satisfy 2:1 ratio"
312 << " when refining cell " << iter.key()
313 << abort(FatalError);
314 }
315 }
316 }
317
318
319 hexRefiner.setRefinement(consistentCells, meshMod);
320
321 // Change mesh, no inflation
322 autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);
323 const mapPolyMesh& morphMap = morphMapPtr();
324
325 if (morphMap.hasMotionPoints())
326 {
327 mesh.movePoints(morphMap.preMotionPoints());
328 }
329
330 if (writeMesh)
331 {
332 mesh.write();
333 }
334
335 if (debug)
336 {
337 Pout<< "multiDirRefinement : updated mesh at time "
338 << mesh.time().timeName() << endl;
339 }
340
341 hexRefiner.updateMesh(morphMap);
342
343 // Collect all cells originating from same old cell (original + 7 extra)
344
345 forAll(consistentCells, i)
346 {
347 addedCells_[consistentCells[i]].setSize(8);
348 }
349 labelList nAddedCells(addedCells_.size(), Zero);
350
351 const labelList& cellMap = morphMap.cellMap();
352
353 forAll(cellMap, celli)
354 {
355 const label oldCelli = cellMap[celli];
356
357 if (addedCells_[oldCelli].size())
358 {
359 addedCells_[oldCelli][nAddedCells[oldCelli]++] = celli;
360 }
361 }
362}
363
364
365void Foam::multiDirRefinement::refineAllDirs
366(
367 polyMesh& mesh,
368 List<vectorField>& cellDirections,
369 const cellLooper& cellWalker,
370 undoableMeshCutter& cutter,
371 const bool writeMesh
372)
373{
374 // Iterator
375 refinementIterator refiner(mesh, cutter, cellWalker, writeMesh);
376
377 forAll(cellDirections, dirI)
378 {
379 if (debug)
380 {
381 Pout<< "multiDirRefinement : Refining " << cellLabels_.size()
382 << " cells in direction " << dirI << endl
383 << endl;
384 }
385
386 const vectorField& dirField = cellDirections[dirI];
387
388 // Combine cell to be cut with direction to cut in.
389 // If dirField is only one element use this for all cells.
390
391 List<refineCell> refCells(cellLabels_.size());
392
393 if (dirField.size() == 1)
394 {
395 // Uniform directions.
396 if (debug)
397 {
398 Pout<< "multiDirRefinement : Uniform refinement:"
399 << dirField[0] << endl;
400 }
401
402 forAll(refCells, refI)
403 {
404 const label celli = cellLabels_[refI];
405
406 refCells[refI] = refineCell(celli, dirField[0]);
407 }
408 }
409 else
410 {
411 // Non uniform directions.
412 forAll(refCells, refI)
413 {
414 const label celli = cellLabels_[refI];
415
416 refCells[refI] = refineCell(celli, dirField[celli]);
417 }
418 }
419
420 // Do refine mesh (multiple iterations). Remember added cells.
421 Map<label> splitMap = refiner.setRefinement(refCells);
422
423 // Update overall map of added cells
424 addCells(mesh, splitMap);
425
426 // Add added cells to list of cells to refine in next iteration
427 addCells(splitMap, cellLabels_);
428
429 // Update refinement direction for added cells.
430 if (dirField.size() != 1)
431 {
432 forAll(cellDirections, i)
433 {
434 update(splitMap, cellDirections[i]);
435 }
436 }
437
438 if (debug)
439 {
440 Pout<< "multiDirRefinement : Done refining direction " << dirI
441 << " resulting in " << cellLabels_.size() << " cells" << nl
442 << endl;
443 }
444 }
445}
446
447
448void Foam::multiDirRefinement::refineFromDict
449(
450 polyMesh& mesh,
451 List<vectorField>& cellDirections,
452 const dictionary& dict,
453 const bool writeMesh
454)
455{
456 // How to walk cell circumference.
457 const bool pureGeomCut(dict.get<bool>("geometricCut"));
458
459 autoPtr<cellLooper> cellWalker;
460 if (pureGeomCut)
461 {
462 cellWalker.reset(new geomCellLooper(mesh));
463 }
464 else
465 {
466 cellWalker.reset(new hexCellLooper(mesh));
467 }
468
469 // Construct undoable refinement topology modifier.
470 //Note: undoability switched off.
471 // Might want to reconsider if needs to be possible. But then can always
472 // use other constructor.
473 undoableMeshCutter cutter(mesh, false);
474
475 refineAllDirs(mesh, cellDirections, cellWalker(), cutter, writeMesh);
476}
477
478
479
480// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
481
482// Construct from dictionary
484(
485 polyMesh& mesh,
486 const labelList& cellLabels, // list of cells to refine
487 const dictionary& dict
488)
489:
490 cellLabels_(cellLabels),
491 addedCells_(mesh.nCells())
492{
493 const bool useHex(dict.get<bool>("useHexTopology"));
494
495 const bool writeMesh(dict.get<bool>("writeMesh"));
496
497 const wordList dirNames(dict.get<wordList>("directions"));
498
499 if (useHex && dirNames.size() == 3)
500 {
501 // Filter out hexes from cellLabels_
502 labelList hexCells(splitOffHex(mesh));
503
504 refineHex8(mesh, hexCells, writeMesh);
505 }
506
507 label nRemainingCells = cellLabels_.size();
508
509 reduce(nRemainingCells, sumOp<label>());
510
511 if (nRemainingCells > 0)
512 {
513 // Any cells to refine using meshCutter
514
515 // Determine directions for every cell. Can either be uniform
516 // (size = 1) or non-uniform (one vector per cell)
517 directions cellDirections(mesh, dict);
518
519 refineFromDict(mesh, cellDirections, dict, writeMesh);
520 }
521}
522
523
524// Construct from dictionary and directions to refine.
526(
527 polyMesh& mesh,
528 const labelList& cellLabels, // list of cells to refine
529 const List<vectorField>& cellDirs, // Explicitly provided directions
530 const dictionary& dict
531)
532:
533 cellLabels_(cellLabels),
534 addedCells_(mesh.nCells())
535{
536 const bool useHex(dict.get<bool>("useHexTopology"));
537
538 const bool writeMesh(dict.get<bool>("writeMesh"));
539
540 const wordList dirNames(dict.get<wordList>("directions"));
541
542 if (useHex && dirNames.size() == 3)
543 {
544 // Filter out hexes from cellLabels_
545 labelList hexCells(splitOffHex(mesh));
546
547 refineHex8(mesh, hexCells, writeMesh);
548 }
549
550 label nRemainingCells = cellLabels_.size();
551
552 reduce(nRemainingCells, sumOp<label>());
553
554 if (nRemainingCells > 0)
555 {
556 // Any cells to refine using meshCutter
557
558 // Working copy of cells to refine
559 List<vectorField> cellDirections(cellDirs);
560
561 refineFromDict(mesh, cellDirections, dict, writeMesh);
562 }
563}
564
565
566// Construct from components. Implies meshCutter since directions provided.
568(
569 polyMesh& mesh,
570 undoableMeshCutter& cutter, // actual mesh modifier
571 const cellLooper& cellWalker, // how to cut a single cell with
572 // a plane
573 const labelList& cellLabels, // list of cells to refine
574 const List<vectorField>& cellDirs,
575 const bool writeMesh // write intermediate meshes
576)
577:
578 cellLabels_(cellLabels),
579 addedCells_(mesh.nCells())
580{
581 // Working copy of cells to refine
582 List<vectorField> cellDirections(cellDirs);
583
584 refineAllDirs(mesh, cellDirections, cellWalker, cutter, writeMesh);
585}
586
587
588// ************************************************************************* //
void setSize(const label n)
Alias for resize()
Definition: List.H:218
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Definition: cellLooper.H:75
reference ref() const
A reference to the entry (Error if not found)
Definition: dictionary.H:225
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Set of directions for each cell in the mesh. Either uniform and size=1 or one set of directions per c...
Definition: directions.H:70
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1079
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
virtual void movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:895
virtual bool update()
Update the mesh for both mesh motion and topology change.
Does multiple pass refinement to refine cells in multiple directions.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:324
const cellShapeList & cellShapes() const
Return cell shapes.
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
The main refinement handler. Gets cellCuts which is structure that describes which cells are to be cu...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
rDeltaTY field()
mesh update()
const cellModel & hex
cellShapeList cellShapes
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Namespace for OpenFOAM.
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:45
List< label > labelList
A List of labels.
Definition: List.H:66
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278