FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/Time.H>
33 #include <meshTools/topoSet.H>
34 #include <dynamicMesh/directions.H>
35 #include <dynamicMesh/hexRef8.H>
36 #include <OpenFOAM/mapPolyMesh.H>
38 #include <OpenFOAM/ListOps.H>
39 #include <OpenFOAM/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  for
97  (
98  Map<label>::const_iterator iter = splitMap.begin();
99  iter != splitMap.end();
100  ++iter
101  )
102  {
103  field[iter()] = field[iter.key()];
104  }
105 }
106 
107 
108 // Add added cells to labelList
109 void Foam::multiDirRefinement::addCells
110 (
111  const Map<label>& splitMap,
112  labelList& labels
113 )
114 {
115  label newCellI = labels.size();
116 
117  labels.setSize(labels.size() + splitMap.size());
118 
119  for
120  (
121  Map<label>::const_iterator iter = splitMap.begin();
122  iter != splitMap.end();
123  ++iter
124  )
125  {
126  labels[newCellI++] = iter();
127  }
128 }
129 
130 
131 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
132 
133 void Foam::multiDirRefinement::addCells
134 (
135  const primitiveMesh& mesh,
136  const Map<label>& splitMap
137 )
138 {
139  // Construct inverse addressing: from new to original cell.
140  labelList origCell(mesh.nCells(), -1);
141 
142  forAll(addedCells_, cellI)
143  {
144  const labelList& added = addedCells_[cellI];
145 
146  forAll(added, i)
147  {
148  label slave = added[i];
149 
150  if (origCell[slave] == -1)
151  {
152  origCell[slave] = cellI;
153  }
154  else if (origCell[slave] != cellI)
155  {
157  (
158  "multiDirRefinement::addCells(const primitiveMesh&"
159  ", const Map<label>&"
160  ) << "Added cell " << slave << " has two different masters:"
161  << origCell[slave] << " , " << cellI
162  << abort(FatalError);
163  }
164  }
165  }
166 
167 
168  for
169  (
170  Map<label>::const_iterator iter = splitMap.begin();
171  iter != splitMap.end();
172  ++iter
173  )
174  {
175  label masterI = iter.key();
176  label newCellI = iter();
177 
178  while (origCell[masterI] != -1 && origCell[masterI] != masterI)
179  {
180  masterI = origCell[masterI];
181  }
182 
183  if (masterI >= addedCells_.size())
184  {
186  (
187  "multiDirRefinement::addCells(const primitiveMesh&"
188  ", const Map<label>&"
189  ) << "Map of added cells contains master cell " << masterI
190  << " which is not a valid cell number" << endl
191  << "This means that the mesh is not consistent with the"
192  << " done refinement" << endl
193  << "newCell:" << newCellI << abort(FatalError);
194  }
195 
196  labelList& added = addedCells_[masterI];
197 
198  if (added.empty())
199  {
200  added.setSize(2);
201  added[0] = masterI;
202  added[1] = newCellI;
203  }
204  else if (findIndex(added, newCellI) == -1)
205  {
206  label sz = added.size();
207  added.setSize(sz + 1);
208  added[sz] = newCellI;
209  }
210  }
211 }
212 
213 
214 Foam::labelList Foam::multiDirRefinement::splitOffHex(const primitiveMesh& mesh)
215 {
216  const cellModel& hex = *(cellModeller::lookup("hex"));
217 
218  const cellShapeList& cellShapes = mesh.cellShapes();
219 
220  // Split cellLabels_ into two lists.
221 
222  labelList nonHexLabels(cellLabels_.size());
223  label nonHexI = 0;
224 
225  labelList hexLabels(cellLabels_.size());
226  label hexI = 0;
227 
228  forAll(cellLabels_, i)
229  {
230  label cellI = cellLabels_[i];
231 
232  if (cellShapes[cellI].model() == hex)
233  {
234  hexLabels[hexI++] = cellI;
235  }
236  else
237  {
238  nonHexLabels[nonHexI++] = cellI;
239  }
240  }
241 
242  nonHexLabels.setSize(nonHexI);
243 
244  cellLabels_.transfer(nonHexLabels);
245 
246  hexLabels.setSize(hexI);
247 
248  return hexLabels;
249 }
250 
251 
252 void Foam::multiDirRefinement::refineHex8
253 (
254  polyMesh& mesh,
255  const labelList& hexCells,
256  const bool writeMesh
257 )
258 {
259  if (debug)
260  {
261  Pout<< "multiDirRefinement : Refining hexes " << hexCells.size()
262  << endl;
263  }
264 
265  hexRef8 hexRefiner
266  (
267  mesh,
268  labelList(mesh.nCells(), 0), // cellLevel
269  labelList(mesh.nPoints(), 0), // pointLevel
270  refinementHistory
271  (
272  IOobject
273  (
274  "refinementHistory",
275  mesh.facesInstance(),
277  mesh,
280  false
281  ),
282  List<refinementHistory::splitCell8>(0),
283  labelList(0)
284  ) // refinement history
285  );
286 
287  polyTopoChange meshMod(mesh);
288 
289  labelList consistentCells
290  (
291  hexRefiner.consistentRefinement
292  (
293  hexCells,
294  true // buffer layer
295  )
296  );
297 
298  // Check that consistentRefinement hasn't added cells
299  {
300  // Create count 1 for original cells
301  Map<label> hexCellSet(2*hexCells.size());
302  forAll(hexCells, i)
303  {
304  hexCellSet.insert(hexCells[i], 1);
305  }
306 
307  // Increment count
308  forAll(consistentCells, i)
309  {
310  label cellI = consistentCells[i];
311 
312  Map<label>::iterator iter = hexCellSet.find(cellI);
313 
314  if (iter == hexCellSet.end())
315  {
317  (
318  "multiDirRefinement::refineHex8"
319  "(polyMesh&, const labelList&, const bool)"
320  ) << "Resulting mesh would not satisfy 2:1 ratio"
321  << " when refining cell " << cellI << abort(FatalError);
322  }
323  else
324  {
325  iter() = 2;
326  }
327  }
328 
329  // Check if all been visited (should always be since
330  // consistentRefinement set up to extend set.
331  forAllConstIter(Map<label>, hexCellSet, iter)
332  {
333  if (iter() != 2)
334  {
336  (
337  "multiDirRefinement::refineHex8"
338  "(polyMesh&, const labelList&, const bool)"
339  ) << "Resulting mesh would not satisfy 2:1 ratio"
340  << " when refining cell " << iter.key()
341  << abort(FatalError);
342  }
343  }
344  }
345 
346 
347  hexRefiner.setRefinement(consistentCells, meshMod);
348 
349  // Change mesh, no inflation
350  autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);
351  const mapPolyMesh& morphMap = morphMapPtr();
352 
353  if (morphMap.hasMotionPoints())
354  {
355  mesh.movePoints(morphMap.preMotionPoints());
356  }
357 
358  if (writeMesh)
359  {
360  mesh.write();
361  }
362 
363  if (debug)
364  {
365  Pout<< "multiDirRefinement : updated mesh at time "
366  << mesh.time().timeName() << endl;
367  }
368 
369  hexRefiner.updateMesh(morphMap);
370 
371  // Collect all cells originating from same old cell (original + 7 extra)
372 
373  forAll(consistentCells, i)
374  {
375  addedCells_[consistentCells[i]].setSize(8);
376  }
377  labelList nAddedCells(addedCells_.size(), 0);
378 
379  const labelList& cellMap = morphMap.cellMap();
380 
381  forAll(cellMap, cellI)
382  {
383  label oldCellI = cellMap[cellI];
384 
385  if (addedCells_[oldCellI].size())
386  {
387  addedCells_[oldCellI][nAddedCells[oldCellI]++] = cellI;
388  }
389  }
390 }
391 
392 
393 void Foam::multiDirRefinement::refineAllDirs
394 (
395  polyMesh& mesh,
396  List<vectorField>& cellDirections,
397  const cellLooper& cellWalker,
398  undoableMeshCutter& cutter,
399  const bool writeMesh
400 )
401 {
402  // Iterator
403  refinementIterator refiner(mesh, cutter, cellWalker, writeMesh);
404 
405  forAll(cellDirections, dirI)
406  {
407  if (debug)
408  {
409  Pout<< "multiDirRefinement : Refining " << cellLabels_.size()
410  << " cells in direction " << dirI << endl
411  << endl;
412  }
413 
414  const vectorField& dirField = cellDirections[dirI];
415 
416  // Combine cell to be cut with direction to cut in.
417  // If dirField is only one element use this for all cells.
418 
419  List<refineCell> refCells(cellLabels_.size());
420 
421  if (dirField.size() == 1)
422  {
423  // Uniform directions.
424  if (debug)
425  {
426  Pout<< "multiDirRefinement : Uniform refinement:"
427  << dirField[0] << endl;
428  }
429 
430  forAll(refCells, refI)
431  {
432  label cellI = cellLabels_[refI];
433 
434  refCells[refI] = refineCell(cellI, dirField[0]);
435  }
436  }
437  else
438  {
439  // Non uniform directions.
440  forAll(refCells, refI)
441  {
442  label cellI = cellLabels_[refI];
443 
444  refCells[refI] = refineCell(cellI, dirField[cellI]);
445  }
446  }
447 
448  // Do refine mesh (multiple iterations). Remember added cells.
449  Map<label> splitMap = refiner.setRefinement(refCells);
450 
451  // Update overall map of added cells
452  addCells(mesh, splitMap);
453 
454  // Add added cells to list of cells to refine in next iteration
455  addCells(splitMap, cellLabels_);
456 
457  // Update refinement direction for added cells.
458  if (dirField.size() != 1)
459  {
460  forAll(cellDirections, i)
461  {
462  update(splitMap, cellDirections[i]);
463  }
464  }
465 
466  if (debug)
467  {
468  Pout<< "multiDirRefinement : Done refining direction " << dirI
469  << " resulting in " << cellLabels_.size() << " cells" << nl
470  << endl;
471  }
472  }
473 }
474 
475 
476 void Foam::multiDirRefinement::refineFromDict
477 (
478  polyMesh& mesh,
479  List<vectorField>& cellDirections,
480  const dictionary& dict,
481  const bool writeMesh
482 )
483 {
484  // How to walk cell circumference.
485  Switch pureGeomCut(dict.lookup("geometricCut"));
486 
487  autoPtr<cellLooper> cellWalker(NULL);
488  if (pureGeomCut)
489  {
490  cellWalker.reset(new geomCellLooper(mesh));
491  }
492  else
493  {
494  cellWalker.reset(new hexCellLooper(mesh));
495  }
496 
497  // Construct undoable refinement topology modifier.
498  //Note: undoability switched off.
499  // Might want to reconsider if needs to be possible. But then can always
500  // use other constructor.
501  undoableMeshCutter cutter(mesh, false);
502 
503  refineAllDirs(mesh, cellDirections, cellWalker(), cutter, writeMesh);
504 }
505 
506 
507 
508 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
509 
510 // Construct from dictionary
511 Foam::multiDirRefinement::multiDirRefinement
512 (
513  polyMesh& mesh,
514  const labelList& cellLabels, // list of cells to refine
515  const dictionary& dict
516 )
517 :
518  cellLabels_(cellLabels),
519  addedCells_(mesh.nCells())
520 {
521  Switch useHex(dict.lookup("useHexTopology"));
522 
523  Switch writeMesh(dict.lookup("writeMesh"));
524 
525  wordList dirNames(dict.lookup("directions"));
526 
527  if (useHex && dirNames.size() == 3)
528  {
529  // Filter out hexes from cellLabels_
530  labelList hexCells(splitOffHex(mesh));
531 
532  refineHex8(mesh, hexCells, writeMesh);
533  }
534 
535  label nRemainingCells = cellLabels_.size();
536 
537  reduce(nRemainingCells, sumOp<label>());
538 
539  if (nRemainingCells > 0)
540  {
541  // Any cells to refine using meshCutter
542 
543  // Determine directions for every cell. Can either be uniform
544  // (size = 1) or non-uniform (one vector per cell)
545  directions cellDirections(mesh, dict);
546 
547  refineFromDict(mesh, cellDirections, dict, writeMesh);
548  }
549 }
550 
551 
552 // Construct from directionary and directions to refine.
553 Foam::multiDirRefinement::multiDirRefinement
554 (
555  polyMesh& mesh,
556  const labelList& cellLabels, // list of cells to refine
557  const List<vectorField>& cellDirs, // Explicitly provided directions
558  const dictionary& dict
559 )
560 :
561  cellLabels_(cellLabels),
562  addedCells_(mesh.nCells())
563 {
564  Switch useHex(dict.lookup("useHexTopology"));
565 
566  Switch writeMesh(dict.lookup("writeMesh"));
567 
568  wordList dirNames(dict.lookup("directions"));
569 
570  if (useHex && dirNames.size() == 3)
571  {
572  // Filter out hexes from cellLabels_
573  labelList hexCells(splitOffHex(mesh));
574 
575  refineHex8(mesh, hexCells, writeMesh);
576  }
577 
578  label nRemainingCells = cellLabels_.size();
579 
580  reduce(nRemainingCells, sumOp<label>());
581 
582  if (nRemainingCells > 0)
583  {
584  // Any cells to refine using meshCutter
585 
586  // Working copy of cells to refine
587  List<vectorField> cellDirections(cellDirs);
588 
589  refineFromDict(mesh, cellDirections, dict, writeMesh);
590  }
591 }
592 
593 
594 // Construct from components. Implies meshCutter since directions provided.
595 Foam::multiDirRefinement::multiDirRefinement
596 (
597  polyMesh& mesh,
598  undoableMeshCutter& cutter, // actual mesh modifier
599  const cellLooper& cellWalker, // how to cut a single cell with
600  // a plane
601  const labelList& cellLabels, // list of cells to refine
602  const List<vectorField>& cellDirs,
603  const bool writeMesh // write intermediate meshes
604 )
605 :
606  cellLabels_(cellLabels),
607  addedCells_(mesh.nCells())
608 {
609  // Working copy of cells to refine
610  List<vectorField> cellDirections(cellDirs);
611 
612  refineAllDirs(mesh, cellDirections, cellWalker, cutter, writeMesh);
613 }
614 
615 
616 // ************************ vim: set sw=4 sts=4 et: ************************ //