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