1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | cfMesh: A library for mesh generation
5 \\ / A nd | Author: Franjo Juretic (franjo.juretic@c-fields.com)
6 \\/ M anipulation | Copyright (C) Creative Fields, Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of cfMesh.
11 cfMesh 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 3 of the License, or (at your
14 option) any later version.
16 cfMesh 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 cfMesh. If not, see <http://www.gnu.org/licenses/>.
26 \*---------------------------------------------------------------------------*/
28 #include "meshOctreeCreator.H"
29 #include "meshOctreeModifier.H"
30 #include "IOdictionary.H"
33 #include "meshOctreeAutomaticRefinement.H"
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 // Private member functions
49 void meshOctreeCreator::setRootCubeSizeAndRefParameters()
51 meshOctreeModifier octreeModifier(octree_);
52 if( octreeModifier.isRootInitialisedAccess() )
56 WarningIn("void meshOctreeCreator::setRootCubeSizeAndRefParameters()")
57 << "meshDict is not available" << endl;
64 readScalar(meshDictPtr_->lookup("maxCellSize"))
67 octreeModifier.searchRangeAccess() = 0.25 * maxSize;
69 const triSurf& surface = octree_.surface();
70 boundBox& rootBox = octreeModifier.rootBoxAccess();
71 const point c = (rootBox.max() + rootBox.min()) / 2.0;
73 scalar size = rootBox.max().x() - rootBox.min().x();
74 if( (rootBox.max().y() - rootBox.min().y()) > size )
75 size = rootBox.max().y() - rootBox.min().y();
76 if( (rootBox.max().z() - rootBox.min().z()) > size )
77 size = rootBox.max().z() - rootBox.min().z();
79 size += 0.5 * maxSize;
87 const scalar lSize = size / Foam::pow(2, label(globalRefLevel_));
89 if( lSize < (maxSize * (1.0-SMALL)) )
92 0.5 * maxSize * Foam::pow(2, label(globalRefLevel_));
93 rootBox.max() = c + point(bbSize, bbSize, bbSize);
94 rootBox.min() = c - point(bbSize, bbSize, bbSize);
101 } while( !finished );
104 if( Pstream::parRun() )
106 label l = globalRefLevel_;
107 reduce(l, maxOp<label>());
111 Info << "Root box " << rootBox << endl;
112 Info << "Requested cell size corresponds to octree level "
113 << label(globalRefLevel_) << endl;
115 octreeModifier.isRootInitialisedAccess() = true;
117 surfRefLevel_ = globalRefLevel_;
119 //- set other refinement parameters
120 //- set boundary ref level
121 if( meshDictPtr_->found("boundaryCellSize") )
123 direction boundaryRefLevel_ = globalRefLevel_;
124 scalar cs(readScalar(meshDictPtr_->lookup("boundaryCellSize")));
125 cs *= (scalingFactor_ + SMALL);
127 octreeModifier.searchRangeAccess() = cs;
134 const scalar lSize = maxSize / Foam::pow(2, addLevel);
144 } while( !finished );
146 boundaryRefLevel_ += addLevel;
149 if( Pstream::parRun() )
151 label l = boundaryRefLevel_;
152 reduce(l, maxOp<label>());
153 boundaryRefLevel_ = l;
156 Info << "Requested boundary cell size corresponds to octree level "
157 << label(boundaryRefLevel_) << endl;
159 surfRefLevel_ = boundaryRefLevel_;
162 //- set patch-wise ref levels
163 if( meshDictPtr_->found("patchCellSize") )
165 patchRefinementList refPatches;
167 if( meshDictPtr_->isDict("patchCellSize") )
169 const dictionary& dict = meshDictPtr_->subDict("patchCellSize");
170 const wordList patchNames = dict.toc();
172 refPatches.setSize(patchNames.size());
175 forAll(patchNames, patchI)
177 if( !dict.isDict(patchNames[patchI]) )
180 const dictionary& patchDict = dict.subDict(patchNames[patchI]);
181 const scalar cs = readScalar(patchDict.lookup("cellSize"));
183 refPatches[counter] = patchRefinement(patchNames[patchI], cs);
187 refPatches.setSize(counter);
191 patchRefinementList prl(meshDictPtr_->lookup("patchCellSize"));
192 refPatches.transfer(prl);
195 forAll(refPatches, patchI)
197 const label regionI = refPatches[patchI].patchInSurface(surface);
202 scalar cs = refPatches[patchI].cellSize();
203 cs *= (scalingFactor_ + SMALL);
210 const scalar lSize = maxSize / Foam::pow(2, addLevel);
220 } while( !finished );
222 if( Pstream::parRun() )
223 reduce(addLevel, maxOp<label>());
225 const direction level = globalRefLevel_ + addLevel;
226 forAll(surface, triI)
229 (surface[triI].region() == regionI) &&
230 (surfRefLevel_[triI] < level)
232 surfRefLevel_[triI] = level;
237 if( meshDictPtr_->found("subsetCellSize") )
239 patchRefinementList refPatches;
241 if( meshDictPtr_->isDict("subsetCellSize") )
243 const dictionary& dict = meshDictPtr_->subDict("subsetCellSize");
244 const wordList patchNames = dict.toc();
246 refPatches.setSize(patchNames.size());
249 forAll(patchNames, patchI)
251 if( !dict.isDict(patchNames[patchI]) )
254 const dictionary& patchDict = dict.subDict(patchNames[patchI]);
255 const scalar cs = readScalar(patchDict.lookup("cellSize"));
257 refPatches[counter] = patchRefinement(patchNames[patchI], cs);
261 refPatches.setSize(counter);
265 patchRefinementList srl(meshDictPtr_->lookup("subsetCellSize"));
266 refPatches.transfer(srl);
269 forAll(refPatches, patchI)
271 const word subsetName = refPatches[patchI].patchName();
273 const label subsetID = surface.facetSubsetIndex(subsetName);
276 Warning << "Surface subset " << subsetName
277 << " does not exist" << endl;
281 scalar cs = refPatches[patchI].cellSize();
282 cs *= (scalingFactor_+ SMALL);
289 const scalar lSize = maxSize / Foam::pow(2, addLevel);
299 } while( !finished );
301 if( Pstream::parRun() )
302 reduce(addLevel, maxOp<label>());
304 labelLongList subsetFaces;
305 surface.facetsInSubset(subsetID, subsetFaces);
306 const direction level = globalRefLevel_ + addLevel;
307 forAll(subsetFaces, tI)
308 if( surfRefLevel_[subsetFaces[tI]] < level )
309 surfRefLevel_[subsetFaces[tI]] = level;
313 if( meshDictPtr_->found("localRefinement") )
315 if( meshDictPtr_->isDict("localRefinement") )
317 const dictionary& dict = meshDictPtr_->subDict("localRefinement");
318 const wordList entries = dict.toc();
320 //- map patch name to its index
321 std::map<word, label> patchToIndex;
322 forAll(surface.patches(), patchI)
323 patchToIndex[surface.patches()[patchI].name()] = patchI;
325 //- map a facet subset name to its index
326 std::map<word, label> setToIndex;
327 DynList<label> setIDs;
328 surface.facetSubsetIndices(setIDs);
330 setToIndex[surface.facetSubsetName(setIDs[i])] = setIDs[i];
332 //- set refinement for these entries
333 forAll(entries, dictI)
335 if( !dict.isDict(entries[dictI]) )
338 const word& pName = entries[dictI];
339 const dictionary& patchDict = dict.subDict(pName);
343 if( patchDict.found("additionalRefinementLevels") )
346 readLabel(patchDict.lookup("additionalRefinementLevels"));
348 else if( patchDict.found("cellSize") )
351 readScalar(patchDict.lookup("cellSize"));
357 const scalar lSize = maxSize / Foam::pow(2, nLevel);
367 } while( !finished );
370 const direction level = globalRefLevel_ + nLevel;
372 if( patchToIndex.find(pName) != patchToIndex.end() )
374 //- patch-based refinement
375 const label patchI = patchToIndex[pName];
377 forAll(surface, triI)
379 if( surface[triI].region() == patchI )
380 surfRefLevel_[triI] =
381 Foam::max(surfRefLevel_[triI], level);
384 if( setToIndex.find(pName) != setToIndex.end() )
386 //- this is a facet subset
387 const label subsetId = setToIndex[pName];
389 labelLongList facetsInSubset;
390 surface.facetsInSubset(subsetId, facetsInSubset);
392 forAll(facetsInSubset, i)
394 const label triI = facetsInSubset[i];
395 surfRefLevel_[triI] =
396 Foam::max(surfRefLevel_[triI], level);
404 void meshOctreeCreator::refineInsideAndUnknownBoxes()
406 const direction cType = meshOctreeCube::INSIDE + meshOctreeCube::UNKNOWN;
408 refineBoxes(globalRefLevel_, cType);
411 void meshOctreeCreator::createOctreeBoxes()
413 //- set root cube size in order to achieve desired maxCellSize
414 Info << "Setting root cube size and refinement parameters" << endl;
415 setRootCubeSizeAndRefParameters();
417 //- refine to required boundary resolution
418 Info << "Refining boundary" << endl;
421 //- perform automatic octree refinement
422 if( !Pstream::parRun() )
424 Info << "Performing automatic refinement" << endl;
425 meshOctreeAutomaticRefinement autoRef(octree_, *meshDictPtr_, false);
428 autoRef.activateHexRefinement();
430 autoRef.automaticRefinement();
433 //- set up inside-outside information
434 createInsideOutsideInformation();
436 //- refine INSIDE and UNKNOWN boxes to to the given cell size
437 refineInsideAndUnknownBoxes();
439 //- refine boxes inside the geometry objects
440 refineBoxesContainedInObjects();
442 //- make sure that INSIDE and UNKNOWN neighbours of DATA boxes
443 //- have the same or higher refinement level
444 refineBoxesNearDataBoxes(1);
446 //- distribute octree such that each processor has the same number
447 //- of leaf boxes which will be used as mesh cells
448 if( Pstream::parRun() )
450 loadDistribution(true);
453 //- delete octree data which is not needed any more
454 // if( Pstream::parRun() )
456 // meshOctreeModifier om(octree_);
457 // om.reduceMemoryConsumption();
461 void meshOctreeCreator::createOctreeWithRefinedBoundary
463 const direction maxLevel,
464 const label nTrianglesInLeaf
467 const triSurf& surface = octree_.surface();
468 surface.facetEdges();
469 surface.edgeFacets();
470 const boundBox& rootBox = octree_.rootBox();
471 meshOctreeModifier octreeModifier(octree_);
472 List<meshOctreeSlot>& slots = octreeModifier.dataSlotsAccess();
476 const LongList<meshOctreeCube*>& leaves =
477 octreeModifier.leavesAccess();
482 # pragma omp parallel reduction(+ : nMarked)
486 meshOctreeSlot* slotPtr = &slots[omp_get_thread_num()];
488 meshOctreeSlot* slotPtr = &slots[0];
492 # pragma omp for schedule(dynamic, 20)
494 forAll(leaves, leafI)
496 meshOctreeCube& oc = *leaves[leafI];
499 octree_.containedTriangles(leafI, ct);
501 if( (oc.level() < maxLevel) && (ct.size() > nTrianglesInLeaf) )
503 oc.refineCube(surface, rootBox, slotPtr);
512 octreeModifier.createListOfLeaves();
517 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
519 } // End namespace Foam
521 // ************************************************************************* //