1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "multiDirRefinement.H"
28 #include <OpenFOAM/polyMesh.H>
29 #include <dynamicMesh/polyTopoChanger.H>
30 #include <OpenFOAM/Time.H>
31 #include <dynamicMesh/undoableMeshCutter.H>
32 #include <dynamicMesh/hexCellLooper.H>
33 #include <dynamicMesh/geomCellLooper.H>
34 #include <meshTools/topoSet.H>
35 #include <dynamicMesh/directions.H>
36 #include <dynamicMesh/hexRef8.H>
37 #include <OpenFOAM/mapPolyMesh.H>
38 #include <dynamicMesh/polyTopoChange.H>
39 #include <OpenFOAM/ListOps.H>
40 #include <OpenFOAM/cellModeller.H>
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 defineTypeNameAndDebug(multiDirRefinement, 0);
50 // * * * * * * * * * * * * * Private Statc Functions * * * * * * * * * * * * //
52 // Update refCells pattern for split cells. Note refCells is current
53 // list of cells to refine (these should all have been refined)
54 void Foam::multiDirRefinement::addCells
56 const Map<label>& splitMap,
57 List<refineCell>& refCells
60 label newRefI = refCells.size();
62 label oldSize = refCells.size();
64 refCells.setSize(newRefI + splitMap.size());
66 for (label refI = 0; refI < oldSize; refI++)
68 const refineCell& refCell = refCells[refI];
70 Map<label>::const_iterator iter = splitMap.find(refCell.cellNo());
72 if (iter == splitMap.end())
76 "multiDirRefinement::addCells(const Map<label>&"
77 ", List<refineCell>&)"
78 ) << "Problem : cannot find added cell for cell "
79 << refCell.cellNo() << abort(FatalError);
82 refCells[newRefI++] = refineCell(iter(), refCell.direction());
87 // Update vectorField for all the new cells. Takes over value of
89 void Foam::multiDirRefinement::update
91 const Map<label>& splitMap,
95 field.setSize(field.size() + splitMap.size());
99 Map<label>::const_iterator iter = splitMap.begin();
100 iter != splitMap.end();
104 field[iter()] = field[iter.key()];
109 // Add added cells to labelList
110 void Foam::multiDirRefinement::addCells
112 const Map<label>& splitMap,
116 label newCellI = labels.size();
118 labels.setSize(labels.size() + splitMap.size());
122 Map<label>::const_iterator iter = splitMap.begin();
123 iter != splitMap.end();
127 labels[newCellI++] = iter();
132 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
134 void Foam::multiDirRefinement::addCells
136 const primitiveMesh& mesh,
137 const Map<label>& splitMap
140 // Construct inverse addressing: from new to original cell.
141 labelList origCell(mesh.nCells(), -1);
143 forAll(addedCells_, cellI)
145 const labelList& added = addedCells_[cellI];
149 label slave = added[i];
151 if (origCell[slave] == -1)
153 origCell[slave] = cellI;
155 else if (origCell[slave] != cellI)
159 "multiDirRefinement::addCells(const primitiveMesh&"
160 ", const Map<label>&"
161 ) << "Added cell " << slave << " has two different masters:"
162 << origCell[slave] << " , " << cellI
163 << abort(FatalError);
171 Map<label>::const_iterator iter = splitMap.begin();
172 iter != splitMap.end();
176 label masterI = iter.key();
177 label newCellI = iter();
179 while (origCell[masterI] != -1 && origCell[masterI] != masterI)
181 masterI = origCell[masterI];
184 if (masterI >= addedCells_.size())
188 "multiDirRefinement::addCells(const primitiveMesh&"
189 ", const Map<label>&"
190 ) << "Map of added cells contains master cell " << masterI
191 << " which is not a valid cell number" << endl
192 << "This means that the mesh is not consistent with the"
193 << " done refinement" << endl
194 << "newCell:" << newCellI << abort(FatalError);
197 labelList& added = addedCells_[masterI];
205 else if (findIndex(added, newCellI) == -1)
207 label sz = added.size();
208 added.setSize(sz + 1);
209 added[sz] = newCellI;
215 Foam::labelList Foam::multiDirRefinement::splitOffHex(const primitiveMesh& mesh)
217 const cellModel& hex = *(cellModeller::lookup("hex"));
219 const cellShapeList& cellShapes = mesh.cellShapes();
221 // Split cellLabels_ into two lists.
223 labelList nonHexLabels(cellLabels_.size());
226 labelList hexLabels(cellLabels_.size());
229 forAll(cellLabels_, i)
231 label cellI = cellLabels_[i];
233 if (cellShapes[cellI].model() == hex)
235 hexLabels[hexI++] = cellI;
239 nonHexLabels[nonHexI++] = cellI;
243 nonHexLabels.setSize(nonHexI);
245 cellLabels_.transfer(nonHexLabels);
247 hexLabels.setSize(hexI);
253 void Foam::multiDirRefinement::refineHex8
256 const labelList& hexCells,
262 Pout<< "multiDirRefinement : Refining hexes " << hexCells.size()
269 labelList(mesh.nCells(), 0), // cellLevel
270 labelList(mesh.nPoints(), 0), // pointLevel
276 mesh.facesInstance(),
277 polyMesh::meshSubDir,
283 List<refinementHistory::splitCell8>(0),
285 ) // refinement history
288 polyTopoChange meshMod(mesh);
290 labelList consistentCells
292 hexRefiner.consistentRefinement
299 // Check that consistentRefinement hasn't added cells
301 // Create count 1 for original cells
302 Map<label> hexCellSet(2*hexCells.size());
305 hexCellSet.insert(hexCells[i], 1);
309 forAll(consistentCells, i)
311 label cellI = consistentCells[i];
313 Map<label>::iterator iter = hexCellSet.find(cellI);
315 if (iter == hexCellSet.end())
319 "multiDirRefinement::refineHex8"
320 "(polyMesh&, const labelList&, const bool)"
321 ) << "Resulting mesh would not satisfy 2:1 ratio"
322 << " when refining cell " << cellI << abort(FatalError);
330 // Check if all been visited (should always be since
331 // consistentRefinement set up to extend set.
332 forAllConstIter(Map<label>, hexCellSet, iter)
338 "multiDirRefinement::refineHex8"
339 "(polyMesh&, const labelList&, const bool)"
340 ) << "Resulting mesh would not satisfy 2:1 ratio"
341 << " when refining cell " << iter.key()
342 << abort(FatalError);
348 hexRefiner.setRefinement(consistentCells, meshMod);
350 // Change mesh, no inflation
351 autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);
352 const mapPolyMesh& morphMap = morphMapPtr();
354 if (morphMap.hasMotionPoints())
356 mesh.movePoints(morphMap.preMotionPoints());
366 Pout<< "multiDirRefinement : updated mesh at time "
367 << mesh.time().timeName() << endl;
370 hexRefiner.updateMesh(morphMap);
372 // Collect all cells originating from same old cell (original + 7 extra)
374 forAll(consistentCells, i)
376 addedCells_[consistentCells[i]].setSize(8);
378 labelList nAddedCells(addedCells_.size(), 0);
380 const labelList& cellMap = morphMap.cellMap();
382 forAll(cellMap, cellI)
384 label oldCellI = cellMap[cellI];
386 if (addedCells_[oldCellI].size())
388 addedCells_[oldCellI][nAddedCells[oldCellI]++] = cellI;
394 void Foam::multiDirRefinement::refineAllDirs
397 List<vectorField>& cellDirections,
398 const cellLooper& cellWalker,
399 undoableMeshCutter& cutter,
404 refinementIterator refiner(mesh, cutter, cellWalker, writeMesh);
406 forAll(cellDirections, dirI)
410 Pout<< "multiDirRefinement : Refining " << cellLabels_.size()
411 << " cells in direction " << dirI << endl
415 const vectorField& dirField = cellDirections[dirI];
417 // Combine cell to be cut with direction to cut in.
418 // If dirField is only one element use this for all cells.
420 List<refineCell> refCells(cellLabels_.size());
422 if (dirField.size() == 1)
424 // Uniform directions.
427 Pout<< "multiDirRefinement : Uniform refinement:"
428 << dirField[0] << endl;
431 forAll(refCells, refI)
433 label cellI = cellLabels_[refI];
435 refCells[refI] = refineCell(cellI, dirField[0]);
440 // Non uniform directions.
441 forAll(refCells, refI)
443 label cellI = cellLabels_[refI];
445 refCells[refI] = refineCell(cellI, dirField[cellI]);
449 // Do refine mesh (multiple iterations). Remember added cells.
450 Map<label> splitMap = refiner.setRefinement(refCells);
452 // Update overall map of added cells
453 addCells(mesh, splitMap);
455 // Add added cells to list of cells to refine in next iteration
456 addCells(splitMap, cellLabels_);
458 // Update refinement direction for added cells.
459 if (dirField.size() != 1)
461 forAll(cellDirections, i)
463 update(splitMap, cellDirections[i]);
469 Pout<< "multiDirRefinement : Done refining direction " << dirI
470 << " resulting in " << cellLabels_.size() << " cells" << nl
477 void Foam::multiDirRefinement::refineFromDict
480 List<vectorField>& cellDirections,
481 const dictionary& dict,
485 // How to walk cell circumference.
486 Switch pureGeomCut(dict.lookup("geometricCut"));
488 autoPtr<cellLooper> cellWalker(NULL);
491 cellWalker.reset(new geomCellLooper(mesh));
495 cellWalker.reset(new hexCellLooper(mesh));
498 // Construct undoable refinement topology modifier.
499 //Note: undoability switched off.
500 // Might want to reconsider if needs to be possible. But then can always
501 // use other constructor.
502 undoableMeshCutter cutter(mesh, false);
504 refineAllDirs(mesh, cellDirections, cellWalker(), cutter, writeMesh);
509 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
511 // Construct from dictionary
512 Foam::multiDirRefinement::multiDirRefinement
515 const labelList& cellLabels, // list of cells to refine
516 const dictionary& dict
519 cellLabels_(cellLabels),
520 addedCells_(mesh.nCells())
522 Switch useHex(dict.lookup("useHexTopology"));
524 Switch writeMesh(dict.lookup("writeMesh"));
526 wordList dirNames(dict.lookup("directions"));
528 if (useHex && dirNames.size() == 3)
530 // Filter out hexes from cellLabels_
531 labelList hexCells(splitOffHex(mesh));
533 refineHex8(mesh, hexCells, writeMesh);
536 label nRemainingCells = cellLabels_.size();
538 reduce(nRemainingCells, sumOp<label>());
540 if (nRemainingCells > 0)
542 // Any cells to refine using meshCutter
544 // Determine directions for every cell. Can either be uniform
545 // (size = 1) or non-uniform (one vector per cell)
546 directions cellDirections(mesh, dict);
548 refineFromDict(mesh, cellDirections, dict, writeMesh);
553 // Construct from directionary and directions to refine.
554 Foam::multiDirRefinement::multiDirRefinement
557 const labelList& cellLabels, // list of cells to refine
558 const List<vectorField>& cellDirs, // Explicitly provided directions
559 const dictionary& dict
562 cellLabels_(cellLabels),
563 addedCells_(mesh.nCells())
565 Switch useHex(dict.lookup("useHexTopology"));
567 Switch writeMesh(dict.lookup("writeMesh"));
569 wordList dirNames(dict.lookup("directions"));
571 if (useHex && dirNames.size() == 3)
573 // Filter out hexes from cellLabels_
574 labelList hexCells(splitOffHex(mesh));
576 refineHex8(mesh, hexCells, writeMesh);
579 label nRemainingCells = cellLabels_.size();
581 reduce(nRemainingCells, sumOp<label>());
583 if (nRemainingCells > 0)
585 // Any cells to refine using meshCutter
587 // Working copy of cells to refine
588 List<vectorField> cellDirections(cellDirs);
590 refineFromDict(mesh, cellDirections, dict, writeMesh);
595 // Construct from components. Implies meshCutter since directions provided.
596 Foam::multiDirRefinement::multiDirRefinement
599 undoableMeshCutter& cutter, // actual mesh modifier
600 const cellLooper& cellWalker, // how to cut a single cell with
602 const labelList& cellLabels, // list of cells to refine
603 const List<vectorField>& cellDirs,
604 const bool writeMesh // write intermediate meshes
607 cellLabels_(cellLabels),
608 addedCells_(mesh.nCells())
610 // Working copy of cells to refine
611 List<vectorField> cellDirections(cellDirs);
613 refineAllDirs(mesh, cellDirections, cellWalker, cutter, writeMesh);
617 // ************************ vim: set sw=4 sts=4 et: ************************ //