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