BUG: UListIO: byteSize overflowing on really big faceLists
[OpenFOAM-2.0.x.git] / src / mesh / autoMesh / autoHexMesh / autoHexMeshDriver / autoRefineDriver.C
blob1021e164322c848813cf96722168ecb05390cd76
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
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
13     the Free Software Foundation, either version 3 of the License, or
14     (at your 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
19     for more details.
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/>.
24 \*----------------------------------------------------------------------------*/
26 #include "autoRefineDriver.H"
27 #include "meshRefinement.H"
28 #include "fvMesh.H"
29 #include "Time.H"
30 #include "cellSet.H"
31 #include "syncTools.H"
32 #include "refinementParameters.H"
33 #include "refinementSurfaces.H"
34 #include "refinementFeatures.H"
35 #include "shellSurfaces.H"
36 #include "mapDistributePolyMesh.H"
37 #include "unitConversion.H"
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 namespace Foam
44 defineTypeNameAndDebug(autoRefineDriver, 0);
46 } // End namespace Foam
49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
52 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
54 // Construct from components
55 Foam::autoRefineDriver::autoRefineDriver
57     meshRefinement& meshRefiner,
58     decompositionMethod& decomposer,
59     fvMeshDistribute& distributor,
60     const labelList& globalToPatch
63     meshRefiner_(meshRefiner),
64     decomposer_(decomposer),
65     distributor_(distributor),
66     globalToPatch_(globalToPatch)
70 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
72 Foam::label Foam::autoRefineDriver::featureEdgeRefine
74     const refinementParameters& refineParams,
75     const label maxIter,
76     const label minRefine
79     const fvMesh& mesh = meshRefiner_.mesh();
81     label iter = 0;
83     if (meshRefiner_.features().size() && maxIter > 0)
84     {
85         for (; iter < maxIter; iter++)
86         {
87             Info<< nl
88                 << "Feature refinement iteration " << iter << nl
89                 << "------------------------------" << nl
90                 << endl;
92             labelList candidateCells
93             (
94                 meshRefiner_.refineCandidates
95                 (
96                     refineParams.keepPoints()[0],    // For now only use one.
97                     refineParams.curvature(),
99                     true,               // featureRefinement
100                     false,              // internalRefinement
101                     false,              // surfaceRefinement
102                     false,              // curvatureRefinement
103                     refineParams.maxGlobalCells(),
104                     refineParams.maxLocalCells()
105                 )
106             );
107             labelList cellsToRefine
108             (
109                 meshRefiner_.meshCutter().consistentRefinement
110                 (
111                     candidateCells,
112                     true
113                 )
114             );
115             Info<< "Determined cells to refine in = "
116                 << mesh.time().cpuTimeIncrement() << " s" << endl;
120             label nCellsToRefine = cellsToRefine.size();
121             reduce(nCellsToRefine, sumOp<label>());
123             Info<< "Selected for feature refinement : " << nCellsToRefine
124                 << " cells (out of " << mesh.globalData().nTotalCells()
125                 << ')' << endl;
127             if (nCellsToRefine <= minRefine)
128             {
129                 Info<< "Stopping refining since too few cells selected."
130                     << nl << endl;
131                 break;
132             }
135             if (debug > 0)
136             {
137                 const_cast<Time&>(mesh.time())++;
138             }
141             if
142             (
143                 returnReduce
144                 (
145                     (mesh.nCells() >= refineParams.maxLocalCells()),
146                     orOp<bool>()
147                 )
148             )
149             {
150                 meshRefiner_.balanceAndRefine
151                 (
152                     "feature refinement iteration " + name(iter),
153                     decomposer_,
154                     distributor_,
155                     cellsToRefine,
156                     refineParams.maxLoadUnbalance()
157                 );
158             }
159             else
160             {
161                 meshRefiner_.refineAndBalance
162                 (
163                     "feature refinement iteration " + name(iter),
164                     decomposer_,
165                     distributor_,
166                     cellsToRefine,
167                     refineParams.maxLoadUnbalance()
168                 );
169             }
170         }
171     }
172     return iter;
176 Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
178     const refinementParameters& refineParams,
179     const label maxIter
182     const fvMesh& mesh = meshRefiner_.mesh();
184     // Determine the maximum refinement level over all surfaces. This
185     // determines the minumum number of surface refinement iterations.
186     label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
188     label iter;
189     for (iter = 0; iter < maxIter; iter++)
190     {
191         Info<< nl
192             << "Surface refinement iteration " << iter << nl
193             << "------------------------------" << nl
194             << endl;
197         // Determine cells to refine
198         // ~~~~~~~~~~~~~~~~~~~~~~~~~
199         // Only look at surface intersections (minLevel and surface curvature),
200         // do not do internal refinement (refinementShells)
202         labelList candidateCells
203         (
204             meshRefiner_.refineCandidates
205             (
206                 refineParams.keepPoints()[0],
207                 refineParams.curvature(),
209                 false,              // featureRefinement
210                 false,              // internalRefinement
211                 true,               // surfaceRefinement
212                 true,               // curvatureRefinement
213                 refineParams.maxGlobalCells(),
214                 refineParams.maxLocalCells()
215             )
216         );
217         labelList cellsToRefine
218         (
219             meshRefiner_.meshCutter().consistentRefinement
220             (
221                 candidateCells,
222                 true
223             )
224         );
225         Info<< "Determined cells to refine in = "
226             << mesh.time().cpuTimeIncrement() << " s" << endl;
229         label nCellsToRefine = cellsToRefine.size();
230         reduce(nCellsToRefine, sumOp<label>());
232         Info<< "Selected for refinement : " << nCellsToRefine
233             << " cells (out of " << mesh.globalData().nTotalCells()
234             << ')' << endl;
236         // Stop when no cells to refine or have done minimum nessecary
237         // iterations and not enough cells to refine.
238         if
239         (
240             nCellsToRefine == 0
241          || (
242                 iter >= overallMaxLevel
243              && nCellsToRefine <= refineParams.minRefineCells()
244             )
245         )
246         {
247             Info<< "Stopping refining since too few cells selected."
248                 << nl << endl;
249             break;
250         }
253         if (debug)
254         {
255             const_cast<Time&>(mesh.time())++;
256         }
259         if
260         (
261             returnReduce
262             (
263                 (mesh.nCells() >= refineParams.maxLocalCells()),
264                 orOp<bool>()
265             )
266         )
267         {
268             meshRefiner_.balanceAndRefine
269             (
270                 "surface refinement iteration " + name(iter),
271                 decomposer_,
272                 distributor_,
273                 cellsToRefine,
274                 refineParams.maxLoadUnbalance()
275             );
276         }
277         else
278         {
279             meshRefiner_.refineAndBalance
280             (
281                 "surface refinement iteration " + name(iter),
282                 decomposer_,
283                 distributor_,
284                 cellsToRefine,
285                 refineParams.maxLoadUnbalance()
286             );
287         }
288     }
289     return iter;
293 void Foam::autoRefineDriver::removeInsideCells
295     const refinementParameters& refineParams,
296     const label nBufferLayers
299     Info<< nl
300         << "Removing mesh beyond surface intersections" << nl
301         << "------------------------------------------" << nl
302         << endl;
304     const fvMesh& mesh = meshRefiner_.mesh();
306     if (debug)
307     {
308        const_cast<Time&>(mesh.time())++;
309     }
311     meshRefiner_.splitMesh
312     (
313         nBufferLayers,                  // nBufferLayers
314         globalToPatch_,
315         refineParams.keepPoints()[0]
316     );
318     if (debug)
319     {
320         Pout<< "Writing subsetted mesh to time "
321             << meshRefiner_.timeName() << '.' << endl;
322         meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
323         Pout<< "Dumped mesh in = "
324             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
325     }
329 Foam::label Foam::autoRefineDriver::shellRefine
331     const refinementParameters& refineParams,
332     const label maxIter
335     const fvMesh& mesh = meshRefiner_.mesh();
337     // Mark current boundary faces with 0. Have meshRefiner maintain them.
338     meshRefiner_.userFaceData().setSize(1);
340     // mark list to remove any refined faces
341     meshRefiner_.userFaceData()[0].first() = meshRefinement::REMOVE;
342     meshRefiner_.userFaceData()[0].second() = createWithValues<labelList>
343     (
344         mesh.nFaces(),
345         -1,
346         meshRefiner_.intersectedFaces(),
347         0
348     );
350     // Determine the maximum refinement level over all volume refinement
351     // regions. This determines the minumum number of shell refinement
352     // iterations.
353     label overallMaxShellLevel = meshRefiner_.shells().maxLevel();
355     label iter;
356     for (iter = 0; iter < maxIter; iter++)
357     {
358         Info<< nl
359             << "Shell refinement iteration " << iter << nl
360             << "----------------------------" << nl
361             << endl;
363         labelList candidateCells
364         (
365             meshRefiner_.refineCandidates
366             (
367                 refineParams.keepPoints()[0],
368                 refineParams.curvature(),
370                 false,              // featureRefinement
371                 true,               // internalRefinement
372                 false,              // surfaceRefinement
373                 false,              // curvatureRefinement
374                 refineParams.maxGlobalCells(),
375                 refineParams.maxLocalCells()
376             )
377         );
379         if (debug)
380         {
381             Pout<< "Dumping " << candidateCells.size()
382                 << " cells to cellSet candidateCellsFromShells." << endl;
384             cellSet c(mesh, "candidateCellsFromShells", candidateCells);
385             c.instance() = meshRefiner_.timeName();
386             c.write();
387         }
389         // Problem choosing starting faces for bufferlayers (bFaces)
390         //  - we can't use the current intersected boundary faces
391         //    (intersectedFaces) since this grows indefinitely
392         //  - if we use 0 faces we don't satisfy bufferLayers from the
393         //    surface.
394         //  - possibly we want to have bFaces only the initial set of faces
395         //    and maintain the list while doing the refinement.
396         labelList bFaces
397         (
398             findIndices(meshRefiner_.userFaceData()[0].second(), 0)
399         );
401         //Info<< "Collected boundary faces : "
402         //    << returnReduce(bFaces.size(), sumOp<label>()) << endl;
404         labelList cellsToRefine;
406         if (refineParams.nBufferLayers() <= 2)
407         {
408             cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement
409             (
410                 refineParams.nBufferLayers(),
411                 candidateCells,                     // cells to refine
412                 bFaces,                             // faces for nBufferLayers
413                 1,                                  // point difference
414                 meshRefiner_.intersectedPoints()    // points to check
415             );
416         }
417         else
418         {
419             cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement2
420             (
421                 refineParams.nBufferLayers(),
422                 candidateCells,                 // cells to refine
423                 bFaces                          // faces for nBufferLayers
424             );
425         }
427         Info<< "Determined cells to refine in = "
428             << mesh.time().cpuTimeIncrement() << " s" << endl;
431         label nCellsToRefine = cellsToRefine.size();
432         reduce(nCellsToRefine, sumOp<label>());
434         Info<< "Selected for internal refinement : " << nCellsToRefine
435             << " cells (out of " << mesh.globalData().nTotalCells()
436             << ')' << endl;
438         // Stop when no cells to refine or have done minimum nessecary
439         // iterations and not enough cells to refine.
440         if
441         (
442             nCellsToRefine == 0
443          || (
444                 iter >= overallMaxShellLevel
445              && nCellsToRefine <= refineParams.minRefineCells()
446             )
447         )
448         {
449             Info<< "Stopping refining since too few cells selected."
450                 << nl << endl;
451             break;
452         }
455         if (debug)
456         {
457             const_cast<Time&>(mesh.time())++;
458         }
460         if
461         (
462             returnReduce
463             (
464                 (mesh.nCells() >= refineParams.maxLocalCells()),
465                 orOp<bool>()
466             )
467         )
468         {
469             meshRefiner_.balanceAndRefine
470             (
471                 "shell refinement iteration " + name(iter),
472                 decomposer_,
473                 distributor_,
474                 cellsToRefine,
475                 refineParams.maxLoadUnbalance()
476             );
477         }
478         else
479         {
480             meshRefiner_.refineAndBalance
481             (
482                 "shell refinement iteration " + name(iter),
483                 decomposer_,
484                 distributor_,
485                 cellsToRefine,
486                 refineParams.maxLoadUnbalance()
487             );
488         }
489     }
490     meshRefiner_.userFaceData().clear();
492     return iter;
496 void Foam::autoRefineDriver::baffleAndSplitMesh
498     const refinementParameters& refineParams,
499     const bool handleSnapProblems,
500     const dictionary& motionDict
503     Info<< nl
504         << "Splitting mesh at surface intersections" << nl
505         << "---------------------------------------" << nl
506         << endl;
508     const fvMesh& mesh = meshRefiner_.mesh();
510     // Introduce baffles at surface intersections. Note:
511     // meshRefiment::surfaceIndex() will
512     // be like boundary face from now on so not coupled anymore.
513     meshRefiner_.baffleAndSplitMesh
514     (
515         handleSnapProblems,             // detect&remove potential snap problem
516         false,                          // perpendicular edge connected cells
517         scalarField(0),                 // per region perpendicular angle
518         !handleSnapProblems,            // merge free standing baffles?
519         motionDict,
520         const_cast<Time&>(mesh.time()),
521         globalToPatch_,
522         refineParams.keepPoints()[0]
523     );
527 void Foam::autoRefineDriver::zonify
529     const refinementParameters& refineParams
532     // Mesh is at its finest. Do zoning
533     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
534     // This puts all faces with intersection across a zoneable surface
535     // into that surface's faceZone. All cells inside faceZone get given the
536     // same cellZone.
538     if (meshRefiner_.surfaces().getNamedSurfaces().size())
539     {
540         Info<< nl
541             << "Introducing zones for interfaces" << nl
542             << "--------------------------------" << nl
543             << endl;
545         const fvMesh& mesh = meshRefiner_.mesh();
547         if (debug)
548         {
549             const_cast<Time&>(mesh.time())++;
550         }
552         meshRefiner_.zonify
553         (
554             refineParams.keepPoints()[0],
555             refineParams.allowFreeStandingZoneFaces()
556         );
558         if (debug)
559         {
560             Pout<< "Writing zoned mesh to time "
561                 << meshRefiner_.timeName() << '.' << endl;
562             meshRefiner_.write
563             (
564                 debug,
565                 mesh.time().path()/meshRefiner_.timeName()
566             );
567         }
569         // Check that all faces are synced
570         meshRefinement::checkCoupledFaceZones(mesh);
571     }
575 void Foam::autoRefineDriver::splitAndMergeBaffles
577     const refinementParameters& refineParams,
578     const bool handleSnapProblems,
579     const dictionary& motionDict
582     Info<< nl
583         << "Handling cells with snap problems" << nl
584         << "---------------------------------" << nl
585         << endl;
587     const fvMesh& mesh = meshRefiner_.mesh();
589     // Introduce baffles and split mesh
590     if (debug)
591     {
592         const_cast<Time&>(mesh.time())++;
593     }
595     const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
597     meshRefiner_.baffleAndSplitMesh
598     (
599         handleSnapProblems,
600         handleSnapProblems,                 // remove perp edge connected cells
601         perpAngle,                          // perp angle
602         false,                              // merge free standing baffles?
603         //true,                               // merge free standing baffles?
604         motionDict,
605         const_cast<Time&>(mesh.time()),
606         globalToPatch_,
607         refineParams.keepPoints()[0]
608     );
610     if (debug)
611     {
612         const_cast<Time&>(mesh.time())++;
613     }
615     // Duplicate points on baffles that are on more than one cell
616     // region. This will help snapping pull them to separate surfaces.
617     meshRefiner_.dupNonManifoldPoints();
620     // Merge all baffles that are still remaining after duplicating points.
621     List<labelPair> couples
622     (
623         meshRefiner_.getDuplicateFaces   // get all baffles
624         (
625             identity(mesh.nFaces()-mesh.nInternalFaces())
626           + mesh.nInternalFaces()
627         )
628     );
630     label nCouples = returnReduce(couples.size(), sumOp<label>());
632     Info<< "Detected unsplittable baffles : "
633         << nCouples << endl;
635     if (nCouples > 0)
636     {
637         // Actually merge baffles. Note: not exactly parallellized. Should
638         // convert baffle faces into processor faces if they resulted
639         // from them.
640         meshRefiner_.mergeBaffles(couples);
642         if (debug)
643         {
644             // Debug:test all is still synced across proc patches
645             meshRefiner_.checkData();
646         }
647         Info<< "Merged free-standing baffles in = "
648             << mesh.time().cpuTimeIncrement() << " s." << endl;
649     }
651     if (debug)
652     {
653         Pout<< "Writing handleProblemCells mesh to time "
654             << meshRefiner_.timeName() << '.' << endl;
655         meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
656     }
660 void Foam::autoRefineDriver::mergePatchFaces
662     const refinementParameters& refineParams,
663     const dictionary& motionDict
666     const fvMesh& mesh = meshRefiner_.mesh();
668     Info<< nl
669         << "Merge refined boundary faces" << nl
670         << "----------------------------" << nl
671         << endl;
673     if (debug)
674     {
675         const_cast<Time&>(mesh.time())++;
676     }
678     meshRefiner_.mergePatchFacesUndo
679     (
680         Foam::cos(degToRad(45.0)),
681         Foam::cos(degToRad(45.0)),
682         meshRefiner_.meshedPatches(),
683         motionDict
684     );
686     if (debug)
687     {
688         meshRefiner_.checkData();
689     }
691     meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
693     if (debug)
694     {
695         meshRefiner_.checkData();
696     }
700 void Foam::autoRefineDriver::doRefine
702     const dictionary& refineDict,
703     const refinementParameters& refineParams,
704     const bool prepareForSnapping,
705     const dictionary& motionDict
708     Info<< nl
709         << "Refinement phase" << nl
710         << "----------------" << nl
711         << endl;
713     const fvMesh& mesh = meshRefiner_.mesh();
715     // Check that all the keep points are inside the mesh.
716     refineParams.findCells(mesh);
718     // Refine around feature edges
719     featureEdgeRefine
720     (
721         refineParams,
722         100,    // maxIter
723         0       // min cells to refine
724     );
726     // Refine based on surface
727     surfaceOnlyRefine
728     (
729         refineParams,
730         100     // maxIter
731     );
733     // Remove cells (a certain distance) beyond surface intersections
734     removeInsideCells
735     (
736         refineParams,
737         1       // nBufferLayers
738     );
740     // Internal mesh refinement
741     shellRefine
742     (
743         refineParams,
744         100    // maxIter
745     );
747     // Introduce baffles at surface intersections. Remove sections unreachable
748     // from keepPoint.
749     baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
751     // Mesh is at its finest. Do optional zoning.
752     zonify(refineParams);
754     // Pull baffles apart
755     splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
757     // Do something about cells with refined faces on the boundary
758     if (prepareForSnapping)
759     {
760         mergePatchFaces(refineParams, motionDict);
761     }
764     if (Pstream::parRun())
765     {
766         Info<< nl
767             << "Doing final balancing" << nl
768             << "---------------------" << nl
769             << endl;
771         //if (debug)
772         //{
773         //    const_cast<Time&>(mesh.time())++;
774         //}
776         // Do final balancing. Keep zoned faces on one processor since the
777         // snap phase will convert them to baffles and this only works for
778         // internal faces.
779         meshRefiner_.balance
780         (
781             true,
782             false,
783             scalarField(mesh.nCells(), 1), // dummy weights
784             decomposer_,
785             distributor_
786         );
787     }
791 // ************************************************************************* //