Moving cfMesh into place. Updated contibutors list
[foam-extend-3.2.git] / src / mesh / cfMesh / meshLibrary / utilities / octrees / meshOctree / meshOctreeCreator / meshOctreeCreatorCreateOctreeBoxes.C
blob1037b1b35d1bd82048e68cea3388fcdcf964e224
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | cfMesh: A library for mesh generation
4    \\    /   O peration     |
5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
7 -------------------------------------------------------------------------------
8 License
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
19     for more details.
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/>.
24 Description
26 \*---------------------------------------------------------------------------*/
28 #include "meshOctreeCreator.H"
29 #include "meshOctreeModifier.H"
30 #include "IOdictionary.H"
31 #include "boundBox.H"
32 #include "triSurf.H"
33 #include "meshOctreeAutomaticRefinement.H"
35 # ifdef USE_OMP
36 #include <omp.h>
37 # endif
39 //#define DEBUGOctree
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 namespace Foam
46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 // Private member functions
49 void meshOctreeCreator::setRootCubeSizeAndRefParameters()
51     meshOctreeModifier octreeModifier(octree_);
52     if( octreeModifier.isRootInitialisedAccess() )
53         return;
54     if( !meshDictPtr_ )
55     {
56         WarningIn("void meshOctreeCreator::setRootCubeSizeAndRefParameters()")
57             << "meshDict is not available" << endl;
58         return;
59     }
61     const scalar maxSize
62     (
63         scalingFactor_ *
64         readScalar(meshDictPtr_->lookup("maxCellSize"))
65     );
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;
81     globalRefLevel_ = 0;
82     bool finished;
83     do
84     {
85         finished = false;
87         const scalar lSize = size / Foam::pow(2, label(globalRefLevel_));
89         if( lSize < (maxSize * (1.0-SMALL)) )
90         {
91             const scalar bbSize =
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);
95             finished = true;
96         }
97         else
98         {
99             ++globalRefLevel_;
100         }
101     } while( !finished );
103     //- exchange data
104     if( Pstream::parRun() )
105     {
106         label l = globalRefLevel_;
107         reduce(l, maxOp<label>());
108         globalRefLevel_ = l;
109     }
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") )
122     {
123         direction boundaryRefLevel_ = globalRefLevel_;
124         scalar cs(readScalar(meshDictPtr_->lookup("boundaryCellSize")));
125         cs *= (scalingFactor_ + SMALL);
127         octreeModifier.searchRangeAccess() = cs;
129         label addLevel(0);
130         do
131         {
132             finished = false;
134             const scalar lSize = maxSize / Foam::pow(2, addLevel);
136             if( lSize <= cs )
137             {
138                 finished = true;
139             }
140             else
141             {
142                 ++addLevel;
143             }
144         } while( !finished );
146         boundaryRefLevel_ += addLevel;
148         //- exchange data
149         if( Pstream::parRun() )
150         {
151             label l = boundaryRefLevel_;
152             reduce(l, maxOp<label>());
153             boundaryRefLevel_ = l;
154         }
156         Info << "Requested boundary cell size corresponds to octree level "
157             << label(boundaryRefLevel_) << endl;
159         surfRefLevel_ = boundaryRefLevel_;
160     }
162     //- set patch-wise ref levels
163     if( meshDictPtr_->found("patchCellSize") )
164     {
165         patchRefinementList refPatches;
167         if( meshDictPtr_->isDict("patchCellSize") )
168         {
169             const dictionary& dict = meshDictPtr_->subDict("patchCellSize");
170             const wordList patchNames = dict.toc();
172             refPatches.setSize(patchNames.size());
173             label counter(0);
175             forAll(patchNames, patchI)
176             {
177                 if( !dict.isDict(patchNames[patchI]) )
178                     continue;
180                 const dictionary& patchDict = dict.subDict(patchNames[patchI]);
181                 const scalar cs = readScalar(patchDict.lookup("cellSize"));
183                 refPatches[counter] = patchRefinement(patchNames[patchI], cs);
184                 ++counter;
185             }
187             refPatches.setSize(counter);
188         }
189         else
190         {
191             patchRefinementList prl(meshDictPtr_->lookup("patchCellSize"));
192             refPatches.transfer(prl);
193         }
195         forAll(refPatches, patchI)
196         {
197             const label regionI = refPatches[patchI].patchInSurface(surface);
199             if( regionI < 0 )
200                 continue;
202             scalar cs = refPatches[patchI].cellSize();
203             cs *= (scalingFactor_ + SMALL);
205             label addLevel(0);
206             do
207             {
208                 finished = false;
210                 const scalar lSize = maxSize / Foam::pow(2, addLevel);
212                 if( lSize <= cs )
213                 {
214                     finished = true;
215                 }
216                 else
217                 {
218                     ++addLevel;
219                 }
220             } while( !finished );
222             if( Pstream::parRun() )
223                 reduce(addLevel, maxOp<label>());
225             const direction level = globalRefLevel_ + addLevel;
226             forAll(surface, triI)
227             {
228                 if(
229                     (surface[triI].region() == regionI) &&
230                     (surfRefLevel_[triI] < level)
231                 )
232                     surfRefLevel_[triI] = level;
233             }
234         }
235     }
237     if( meshDictPtr_->found("subsetCellSize") )
238     {
239         patchRefinementList refPatches;
241         if( meshDictPtr_->isDict("subsetCellSize") )
242         {
243             const dictionary& dict = meshDictPtr_->subDict("subsetCellSize");
244             const wordList patchNames = dict.toc();
246             refPatches.setSize(patchNames.size());
247             label counter(0);
249             forAll(patchNames, patchI)
250             {
251                 if( !dict.isDict(patchNames[patchI]) )
252                     continue;
254                 const dictionary& patchDict = dict.subDict(patchNames[patchI]);
255                 const scalar cs = readScalar(patchDict.lookup("cellSize"));
257                 refPatches[counter] = patchRefinement(patchNames[patchI], cs);
258                 ++counter;
259             }
261             refPatches.setSize(counter);
262         }
263         else
264         {
265             patchRefinementList srl(meshDictPtr_->lookup("subsetCellSize"));
266             refPatches.transfer(srl);
267         }
269         forAll(refPatches, patchI)
270         {
271             const word subsetName = refPatches[patchI].patchName();
273             const label subsetID = surface.facetSubsetIndex(subsetName);
274             if( subsetID < 0 )
275             {
276                 Warning << "Surface subset " << subsetName
277                     << " does not exist" << endl;
278                 continue;
279             }
281             scalar cs = refPatches[patchI].cellSize();
282             cs *= (scalingFactor_+ SMALL);
284             label addLevel(0);
285             do
286             {
287                 finished = false;
289                 const scalar lSize = maxSize / Foam::pow(2, addLevel);
291                 if( lSize <= cs )
292                 {
293                     finished = true;
294                 }
295                 else
296                 {
297                     ++addLevel;
298                 }
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;
310         }
311     }
313     if( meshDictPtr_->found("localRefinement") )
314     {
315         if( meshDictPtr_->isDict("localRefinement") )
316         {
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);
329             forAll(setIDs, i)
330                 setToIndex[surface.facetSubsetName(setIDs[i])] = setIDs[i];
332             //- set refinement for these entries
333             forAll(entries, dictI)
334             {
335                 if( !dict.isDict(entries[dictI]) )
336                     continue;
338                 const word& pName = entries[dictI];
339                 const dictionary& patchDict = dict.subDict(pName);
341                 label nLevel(0);
343                 if( patchDict.found("additionalRefinementLevels") )
344                 {
345                     nLevel =
346                         readLabel(patchDict.lookup("additionalRefinementLevels"));
347                 }
348                 else if( patchDict.found("cellSize") )
349                 {
350                     const scalar cs =
351                         readScalar(patchDict.lookup("cellSize"));
353                     do
354                     {
355                         finished = false;
357                         const scalar lSize = maxSize / Foam::pow(2, nLevel);
359                         if( lSize <= cs )
360                         {
361                             finished = true;
362                         }
363                         else
364                         {
365                             ++nLevel;
366                         }
367                     } while( !finished );
368                 }
370                 const direction level = globalRefLevel_ + nLevel;
372                 if( patchToIndex.find(pName) != patchToIndex.end() )
373                 {
374                     //- patch-based refinement
375                     const label patchI = patchToIndex[pName];
377                     forAll(surface, triI)
378                     {
379                         if( surface[triI].region() == patchI )
380                             surfRefLevel_[triI] =
381                                 Foam::max(surfRefLevel_[triI], level);
382                     }
383                 }
384                 if( setToIndex.find(pName) != setToIndex.end() )
385                 {
386                     //- this is a facet subset
387                     const label subsetId = setToIndex[pName];
389                     labelLongList facetsInSubset;
390                     surface.facetsInSubset(subsetId, facetsInSubset);
392                     forAll(facetsInSubset, i)
393                     {
394                         const label triI = facetsInSubset[i];
395                         surfRefLevel_[triI] =
396                             Foam::max(surfRefLevel_[triI], level);
397                     }
398                 }
399             }
400         }
401     }
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;
419     refineBoundary();
421     //- perform automatic octree refinement
422     if( !Pstream::parRun() )
423     {
424         Info << "Performing automatic refinement" << endl;
425         meshOctreeAutomaticRefinement autoRef(octree_, *meshDictPtr_, false);
427         if( hexRefinement_ )
428             autoRef.activateHexRefinement();
430         autoRef.automaticRefinement();
431     }
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() )
449     {
450         loadDistribution(true);
451     }
453     //- delete octree data which is not needed any more
454 //    if( Pstream::parRun() )
455 //    {
456 //        meshOctreeModifier om(octree_);
457 //        om.reduceMemoryConsumption();
458 //    }
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();
474     while( true )
475     {
476         const LongList<meshOctreeCube*>& leaves =
477             octreeModifier.leavesAccess();
479         label nMarked(0);
481         # ifdef USE_OMP
482         # pragma omp parallel reduction(+ : nMarked)
483         # endif
484         {
485             # ifdef USE_OMP
486             meshOctreeSlot* slotPtr = &slots[omp_get_thread_num()];
487             # else
488             meshOctreeSlot* slotPtr = &slots[0];
489             # endif
491             # ifdef USE_OMP
492             # pragma omp for schedule(dynamic, 20)
493             # endif
494             forAll(leaves, leafI)
495             {
496                 meshOctreeCube& oc = *leaves[leafI];
498                 DynList<label> ct;
499                 octree_.containedTriangles(leafI, ct);
501                 if( (oc.level() < maxLevel) && (ct.size() > nTrianglesInLeaf) )
502                 {
503                     oc.refineCube(surface, rootBox, slotPtr);
504                     ++nMarked;
505                 }
506             }
507         }
509         if( nMarked == 0 )
510             break;
512         octreeModifier.createListOfLeaves();
514     }
517 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
519 } // End namespace Foam
521 // ************************************************************************* //