gridcreator.hh 60.5 KB
Newer Older
Timo Koch's avatar
Timo Koch committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup InputOutput
Timo Koch's avatar
Timo Koch committed
22
23
 * \brief Provides a grid creator for all supported grid managers with
 *        input file interfaces.
24
25
 *
 * \todo add Petrel grids with dune-cornerpoint
Timo Koch's avatar
Timo Koch committed
26
27
28
29
30
31
32
33
34
35
 */
#ifndef DUMUX_GRID_CREATOR_HH
#define DUMUX_GRID_CREATOR_HH

#include <array>
#include <bitset>
#include <memory>
#include <sstream>

#include <dune/common/exceptions.hh>
36
#include <dune/common/classname.hh>
Timo Koch's avatar
Timo Koch committed
37
38
39
40
41
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/io/file/dgfparser/dgfparser.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/common/gridfactory.hh>
42
#include <dune/grid/common/datahandleif.hh>
Timo Koch's avatar
Timo Koch committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include <dune/grid/utility/structuredgridfactory.hh>

// YaspGrid specific includes
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>

 // OneDGrid specific includes
#include <dune/grid/onedgrid.hh>
#include <dune/grid/io/file/dgfparser/dgfoned.hh>

// UGGrid specific includes
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#endif

// ALUGrid specific includes
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/dgf.hh>
#endif

// FoamGrid specific includes
#if HAVE_DUNE_FOAMGRID
#include <dune/foamgrid/foamgrid.hh>
68
69
#include <dune/foamgrid/dgffoam.hh>
#endif
Timo Koch's avatar
Timo Koch committed
70

71
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
72
#include <dumux/common/parameters.hh>
73
#include <dumux/discretization/methods.hh>
Timo Koch's avatar
Timo Koch committed
74
75
76
77

namespace Dumux
{

78
template<class GridPtr, class GridCreator>
79
struct GridDataHandle : public Dune::CommDataHandleIF<GridDataHandle<GridPtr, GridCreator>, int>
80
81
{
    using GridType = typename GridPtr::element_type;
82
    GridDataHandle( GridPtr& gridPtr)
83
84
85
86
87
88
89
90
    :  gridPtr_(gridPtr), idSet_(gridPtr->localIdSet())
    {
        auto gridView = gridPtr_->levelGridView(0);

        for( const auto& element : elements(gridView, Dune::Partitions::interior))
           std::swap(GridCreator::elementMarkers_[GridCreator::gridFactory().insertionIndex(element)], elData_[idSet_.id(element)]);
    }

91
    ~GridDataHandle()
92
93
94
95
96
97
98
99
100
101
    {
        auto gridView = gridPtr_->levelGridView(0);
        const auto& indexSet = gridView.indexSet();

        GridCreator::elementMarkers_.resize(indexSet.size(0));

        for(const auto& element : elements(gridView))
           std::swap(GridCreator::elementMarkers_[indexSet.index(element)], elData_[idSet_.id(element)]);
    }

102
    Dune::CommDataHandleIF<GridDataHandle<GridPtr, GridCreator>, int>& interface()
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
    { return *this; }

    bool contains (int dim, int codim) const
    { return codim==0; }

    bool fixedSize (int dim, int codim) const
    { return false; }

    template<class EntityType>
    size_t size (const EntityType& e) const
    { return GridCreator::elementMarkers_.size(); }

    template<class MessageBufferImp, class EntityType>
    void gather (MessageBufferImp& buff, const EntityType& e) const
    { buff.write(elData_[idSet_.id(e)]); }

    template<class MessageBufferImp, class EntityType>
    void scatter (MessageBufferImp& buff, const EntityType& e, size_t n)
    { buff.read(elData_[idSet_.id(e)]); }

private:
    GridPtr &gridPtr_;
    using IdSet = typename GridType::LocalIdSet;
    const IdSet &idSet_;
    mutable std::map< typename IdSet::IdType, int> elData_;
};


Timo Koch's avatar
Timo Koch committed
131
/*!
132
 * \ingroup InputOutput
Timo Koch's avatar
Timo Koch committed
133
134
135
 * \brief Provides the grid creator base interface (public) and methods common
 *        to most grid creator specializations (protected).
 */
136
template <class Grid>
Timo Koch's avatar
Timo Koch committed
137
138
class GridCreatorBase
{
139
140
    using Intersection = typename Grid::LeafIntersection;
    using Element = typename Grid::template Codim<0>::Entity;
141
    friend class GridDataHandle<std::shared_ptr<Grid>, GridCreatorBase>;
Timo Koch's avatar
Timo Koch committed
142
public:
143
144
145
146
147
148
149
150
    /*!
     * \brief Make the grid. Implement this method in the specialization of this class for a grid type.
     */
    static void makeGrid(const std::string& modelParamGroup = "")
    {
        DUNE_THROW(Dune::NotImplemented,
            "The GridCreator for grid type " << Dune::className<Grid>() << " is not implemented! Consider providing your own GridCreator.");
    }
Timo Koch's avatar
Timo Koch committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

    /*!
     * \brief Returns a reference to the grid.
     */
    static Grid &grid()
    {
        if(enableDgfGridPointer_)
            return *dgfGridPtr();
        else
            return *gridPtr();
    }

    /*!
     * \brief Call the parameters function of the DGF grid pointer if available
     */
    template <class Entity>
    static const std::vector<double>& parameters(const Entity& entity)
    {
        if(enableDgfGridPointer_)
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        {
            if (entity.hasFather())
            {
                auto level0entity = entity;
                while(level0entity.hasFather())
                    level0entity = level0entity.father();


                return dgfGridPtr().parameters(level0entity);
            }
            else
            {
                return dgfGridPtr().parameters(entity);
            }
        }
Timo Koch's avatar
Timo Koch committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
        else
            DUNE_THROW(Dune::InvalidStateException, "The parameters method is only available if the grid was constructed with a DGF file!");
    }

    /*!
     * \brief Call the parameters function of the DGF grid pointer if available
     */
    template <class GridImp, class IntersectionImp>
    static const Dune::DGFBoundaryParameter::type& parameters(const Dune::Intersection<GridImp, IntersectionImp>& intersection)
    {
        if(enableDgfGridPointer_)
            return dgfGridPtr().parameters(intersection);
        else
            DUNE_THROW(Dune::InvalidStateException, "The parameters method is only available if the grid was constructed with a DGF file!");
    }

201
    /*!
202
203
204
     * \brief Return the boundary domain marker (Gmsh physical entity number) of an intersection
              Only available when using Gmsh with GridParameterGroup.DomainMarkers = 1.
     * \param boundarySegmentIndex The boundary segment index of the intersection (intersection.boundarySegmentIndex()
205
     */
206
    static int getBoundaryDomainMarker(int boundarySegmentIndex)
207
208
    {
        if(enableGmshDomainMarkers_)
209
210
211
        {
            if (boundarySegmentIndex >= grid().numBoundarySegments())
                DUNE_THROW(Dune::RangeError, "Boundary segment index "<< boundarySegmentIndex << " bigger than number of bonudary segments in grid!");
212
            return boundaryMarkers_[boundarySegmentIndex];
213
        }
214
215
216
        else
            DUNE_THROW(Dune::InvalidStateException, "The getBoundaryDomainMarker method is only available if DomainMarkers for Gmsh were enabled!"
                                                     << " If your Gmsh file contains domain markers / physical entities,"
217
                                                     << " enable them by setting 'Grid.DomainMarkers = true' in the input file.");
218
219
220
    }

    /*!
221
222
223
     * \brief Return the element domain marker (Gmsh physical entity number) of an element.
              Only available when using Gmsh with GridParameterGroup.DomainMarkers = 1.
     * \param elementIdx The element index
224
     */
225
    DUNE_DEPRECATED_MSG("This will often produce wrong parameters in case the grid implementation resorts the elements after insertion. Use getElementDomainMarker(element) instead!")
226
    static int getElementDomainMarker(int elementIdx)
227
228
229
230
231
232
233
234
235
236
    {
        if(enableGmshDomainMarkers_)
        {
            if(elementIdx >= grid().levelGridView(0).size(0))
                DUNE_THROW(Dune::RangeError, "Requested element index is bigger than the number of level 0 elements!");
            return elementMarkers_[elementIdx];
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "The getElementDomainMarker method is only available if DomainMarkers for Gmsh were enabled!"
                                                     << " If your Gmsh file contains domain markers / physical entities,"
237
                                                     << " enable them by setting 'Grid.DomainMarkers = true' in the input file.");
238
239
    }

240
241
242
243
244
245
246
247
248
    /*!
     * \brief Return the element domain marker (Gmsh physical entity number) of an element.
              Only available when using Gmsh with GridParameterGroup.DomainMarkers = 1.
     * \param elementIdx The element index
     */
    static int getElementDomainMarker(const Element& element)
    {
        if(enableGmshDomainMarkers_)
        {
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
            // parameters are only given for level 0 elements
            if (element.hasFather())
            {
                auto level0element = element;
                while(level0element.hasFather())
                    level0element = level0element.father();

                // in the parallel case the data is load balanced and then accessed with indices of the index set
                if (grid().comm().size() > 1)
                    return elementMarkers_[gridPtr()->levelGridView(0).indexSet().index(level0element)];
                else
                    return elementMarkers_[gridFactory().insertionIndex(level0element)];
            }
            else
            {
                // in the parallel case the data is load balanced and then accessed with indices of the index set
                if (grid().comm().size() > 1)
                    return elementMarkers_[gridPtr()->levelGridView(0).indexSet().index(element)];
                else
                    return elementMarkers_[gridFactory().insertionIndex(element)];
            }
270
271
272
273
274
275
276
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "The getElementDomainMarker method is only available if DomainMarkers for Gmsh were enabled!"
                                                    << " If your Gmsh file contains domain markers / physical entities,"
                                                    << " enable them by setting 'Grid.DomainMarkers = true' in the input file.");
    }

Timo Koch's avatar
Timo Koch committed
277
278
279
280
281
    /*!
     * \brief Call loadBalance() function of the grid.
     */
    static void loadBalance()
    {
282
        // if we may have dgf parameters use load balancing of the dgf pointer
Timo Koch's avatar
Timo Koch committed
283
        if(enableDgfGridPointer_)
284
        {
Timo Koch's avatar
Timo Koch committed
285
            dgfGridPtr().loadBalance();
286
287
288
289
        }
        // if we have gmsh parameters we have to manually load balance the data
        else if (enableGmshDomainMarkers_)
        {
290
            GridDataHandle<std::shared_ptr<Grid>, GridCreatorBase<Grid>> dh(gridPtr());
291
292
293
            gridPtr()->loadBalance(dh.interface());
            gridPtr()->communicate(dh.interface(), Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
        }
Timo Koch's avatar
Timo Koch committed
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
        else
            gridPtr()->loadBalance();
    }

protected:

    /*!
     * \brief Returns a reference to the grid pointer (std::shared_ptr<Grid>)
     */
    static std::shared_ptr<Grid> &gridPtr()
    {
        if(!enableDgfGridPointer_)
        {
            static std::shared_ptr<Grid> gridPtr_;
            return gridPtr_;
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "You are using DGF. To get the grid pointer use method dgfGridPtr()!");
    }

314
315
316
317
318
319
320
321
322
    /*!
     * \brief Returns a reference to the grid factory
     */
    static Dune::GridFactory<Grid> &gridFactory()
    {
        static Dune::GridFactory<Grid> gridFactory_;
        return gridFactory_;
    }

Timo Koch's avatar
Timo Koch committed
323
    /*!
324
     * \brief Returns a reference to the DGF grid pointer (Dune::GridPtr<Grid>).
Timo Koch's avatar
Timo Koch committed
325
326
327
328
329
330
331
332
333
     */
    static Dune::GridPtr<Grid> &dgfGridPtr()
    {
        if(enableDgfGridPointer_)
        {
            static Dune::GridPtr<Grid> dgfGridPtr_;
            return dgfGridPtr_;
        }
        else
334
            DUNE_THROW(Dune::InvalidStateException, "The DGF grid pointer is only available if the grid was constructed with a DGF file!");
Timo Koch's avatar
Timo Koch committed
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
    }

    /*!
     * \brief Returns the filename extension of a given filename
     */
    static std::string getFileExtension(const std::string& fileName)
    {
        std::size_t i = fileName.rfind('.', fileName.length());
        if (i != std::string::npos)
        {
            return(fileName.substr(i+1, fileName.length() - i));
        }
        else
        {
            DUNE_THROW(Dune::IOError, "Please provide and extension for your grid file ('"<< fileName << "')!");
        }
        return "";
    }

    /*!
     * \brief Makes a grid from a file. We currently support *.dgf (Dune Grid Format) and *.msh (Gmsh mesh format).
     */
357
    static void makeGridFromFile(const std::string& fileName,
358
                                 const std::string& modelParamGroup)
Timo Koch's avatar
Timo Koch committed
359
360
361
362
    {
        // We found a file in the input file...does it have a supported extension?
        const std::string extension = getFileExtension(fileName);
        if(extension != "dgf" && extension != "msh")
363
            DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " only supports DGF (*.dgf) and Gmsh (*.msh) grid files but the specified filename has extension: *."<< extension);
Timo Koch's avatar
Timo Koch committed
364
365
366
367
368
369
370
371
372
373

        // make the grid
        if(extension == "dgf")
        {
            enableDgfGridPointer_ = true;
            dgfGridPtr() = Dune::GridPtr<Grid>(fileName.c_str(), Dune::MPIHelper::getCommunicator());
        }
        if(extension == "msh")
        {
            // get some optional parameters
374
375
376
            const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
            const bool boundarySegments = getParamFromGroup<bool>(modelParamGroup, "Grid.BoundarySegments", false);
            const bool domainMarkers = getParamFromGroup<bool>(modelParamGroup, "Grid.DomainMarkers", false);
377

378
379
380
381
            if (domainMarkers)
                enableGmshDomainMarkers_ = true;

            // as default read it on all processes in parallel
382
383
            if(domainMarkers)
            {
384
385
                Dune::GmshReader<Grid>::read(gridFactory(), fileName, boundaryMarkers_, elementMarkers_, verbose, boundarySegments);
                gridPtr() = std::shared_ptr<Grid>(gridFactory().createGrid());
386
387
388
            }
            else
            {
389
390
                Dune::GmshReader<Grid>::read(gridFactory(), fileName, verbose, boundarySegments);
                gridPtr() = std::shared_ptr<Grid>(gridFactory().createGrid());
391
            }
Timo Koch's avatar
Timo Koch committed
392
393
394
395
396
397
        }
    }

    /*!
     * \brief Makes a grid from a DGF file. This is used by grid managers that only support DGF.
     */
398
    static void makeGridFromDgfFile(const std::string& fileName)
Timo Koch's avatar
Timo Koch committed
399
400
401
402
    {
        // We found a file in the input file...does it have a supported extension?
        const std::string extension = getFileExtension(fileName);
        if(extension != "dgf")
403
            DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " only supports DGF (*.dgf) but the specified filename has extension: *."<< extension);
Timo Koch's avatar
Timo Koch committed
404
405
406
407
408
409
410
411
412
413
414
415
416
417

        enableDgfGridPointer_ = true;
        dgfGridPtr() = Dune::GridPtr<Grid>(fileName.c_str(), Dune::MPIHelper::getCommunicator());
    }

    /*!
     * \brief The cell types for structured grids
     */
    enum CellType {Simplex, Cube};

    /*!
     * \brief Makes a structured cube grid using the structured grid factory
     */
    template <int dim, int dimworld>
418
    static void makeStructuredGrid(CellType cellType,
419
                                   const std::string& modelParamGroup)
Timo Koch's avatar
Timo Koch committed
420
    {
421
422
423
        using GlobalPosition = Dune::FieldVector<typename Grid::ctype, dimworld>;
        const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight");
        const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0));
424

425
426
427
        using CellArray = std::array<unsigned int, dim>;
        CellArray cells; cells.fill(1);
        cells = getParamFromGroup<CellArray>(modelParamGroup, "Grid.Cells", cells);
Timo Koch's avatar
Timo Koch committed
428
429
430
431

        // make the grid
        if (cellType == CellType::Cube)
        {
432
            gridPtr() = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, cells);
Timo Koch's avatar
Timo Koch committed
433
        }
434
        else if (cellType == CellType::Simplex)
Timo Koch's avatar
Timo Koch committed
435
        {
436
            gridPtr() = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, cells);
Timo Koch's avatar
Timo Koch committed
437
        }
438
439
440
441
        else
        {
            DUNE_THROW(Dune::GridError, "Unknown cell type for making structured grid! Choose Cube or Simplex.");
        }
Timo Koch's avatar
Timo Koch committed
442
443
444
445
446
    }

    /*!
     * \brief Refines a grid after construction if GridParameterGroup.Refinement is set in the input file
     */
447
    static void maybeRefineGrid(const std::string& modelParamGroup)
Timo Koch's avatar
Timo Koch committed
448
    {
449
450
        if (haveParamInGroup(modelParamGroup, "Grid.Refinement"))
            grid().globalRefine(getParamFromGroup<int>(modelParamGroup, "Grid.Refinement"));
Timo Koch's avatar
Timo Koch committed
451
452
453
454
455
456
457
    }

    /*!
    * \brief A state variable if the DGF Dune::GridPtr has been enabled.
    *        It is always enabled if a DGF grid file was used to create the grid.
    */
    static bool enableDgfGridPointer_;
458
459

    /*!
460
    * \brief A state variable if domain markers have been read from a Gmsh file.
461
462
463
464
    */
    static bool enableGmshDomainMarkers_;

    /*!
465
466
    * \brief Element and boundary domain markers obtained from Gmsh physical entities
    *        They map from element indices / boundary ids to the physical entity number
467
468
469
    */
    static std::vector<int> elementMarkers_;
    static std::vector<int> boundaryMarkers_;
Timo Koch's avatar
Timo Koch committed
470
471
};

472
473
template <class Grid>
bool GridCreatorBase<Grid>::enableDgfGridPointer_ = false;
Timo Koch's avatar
Timo Koch committed
474

475
476
template <class Grid>
bool GridCreatorBase<Grid>::enableGmshDomainMarkers_ = false;
477

478
479
template <class Grid>
std::vector<int> GridCreatorBase<Grid>::elementMarkers_;
480

481
482
template <class Grid>
std::vector<int> GridCreatorBase<Grid>::boundaryMarkers_;
483

Timo Koch's avatar
Timo Koch committed
484
485
486
487
488
/*!
 * \brief Provides the grid creator implementation for all supported grid managers that constructs a grid
 *        from information in the input file. This class is specialised below for all
 *        supported grid managers. It inherits the functionality of the base class.
 */
489
template <class Grid, DiscretizationMethod discMethod>
490
class GridCreatorImpl : public GridCreatorBase<Grid> {};
Timo Koch's avatar
Timo Koch committed
491
492
493
494
495

/*!
 * \brief Provides the grid creator (this is the class called by the user) for all supported grid managers that constructs a grid
 *        from information in the input file. This class is specialised below for all
 *        supported grid managers. It inherits the functionality of the base class.
496
497
 * \todo  TODO The grid creator is independent of TypeTag now,
 *        it would only need two template parameters, none of the functions use a TypeTag directly
498
 * \todo  This shouldn't depend on FVGridGeometry, think about how to remove discMethod here, too
Timo Koch's avatar
Timo Koch committed
499
500
 */
template <class TypeTag>
501
using GridCreator = GridCreatorImpl<typename GET_PROP_TYPE(TypeTag, Grid), GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
Timo Koch's avatar
Timo Koch committed
502
503
504
505
506

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Specializations //////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////

507
508
/*!
 * \brief Helper class for determining the default overlap in case of parallel yasp grids
509
 * \note the default of 1 works for all overlapping implementation like the cell-centered discretization schemes
510
 */
511
template <DiscretizationMethod discMethod>
512
struct YaspOverlapHelper
513
{
514
    static int getOverlap(const std::string& modelParamGroup)
515
    {
516
517
518
        int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1);
        if (overlap < 1)
            DUNE_THROW(Dune::InvalidStateException, "Parallel non-overlapping grids for cc discretization is not allowed.");
519
520
        return overlap;
    }
521
};
522

523
524
//! specialization for the box method
template <>
525
struct YaspOverlapHelper<DiscretizationMethod::box>
526
527
528
{
    static int getOverlap(const std::string& modelParamGroup)
    { return 0; }
529
530
};

Timo Koch's avatar
Timo Koch committed
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
/*!
 * \brief Provides a grid creator for YaspGrids
 *        from information in the input file
 *
 * All keys are expected to be in group GridParameterGroup.
 * The following keys are recognized:
 * - File : a DGF file to load the coarse grid from
 * - UpperRight : extension of the domain
 * - Cells : the number of cells in each direction
 * - Periodic : true or false for each direction
 * - Overlap : overlap size in cells
 * - Partitioning : a non-standard load-balancing, number of processors per direction
 * - KeepPyhsicalOverlap : whether to keep the physical overlap
 *     in physical size or in number of cells upon refinement
 * - Refinement : the number of global refines to apply initially.
 *
 */
548
549
template<DiscretizationMethod discMethod, class ct, int dim>
class GridCreatorImpl<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> >, discMethod>
550
          : public GridCreatorBase<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> > >
Timo Koch's avatar
Timo Koch committed
551
552
{
public:
553
554
    using Grid = typename Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> >;
    using ParentType = GridCreatorBase<Grid>;
Timo Koch's avatar
Timo Koch committed
555
556
557
558

    /*!
     * \brief Make the grid. This is implemented by specializations of this method.
     */
559
    static void makeGrid(const std::string& modelParamGroup = "")
Timo Koch's avatar
Timo Koch committed
560
561
    {
        // First try to create it from a DGF file in GridParameterGroup.File
562
        if (haveParamInGroup(modelParamGroup, "Grid.File"))
563
        {
564
565
            ParentType::makeGridFromDgfFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"));
            postProcessing_(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
566
567
568
569
            return;
        }

        // Then look for the necessary keys to construct from the input file
570
571
        else if (haveParamInGroup(modelParamGroup, "Grid.UpperRight"))
        {
Timo Koch's avatar
Timo Koch committed
572

573
574
            // get the upper right corner coordinates
            const auto upperRight = getParamFromGroup<Dune::FieldVector<ct, dim>>(modelParamGroup, "Grid.UpperRight");
575

576
577
578
            // number of cells in each direction
            std::array<int, dim> cells; cells.fill(1);
            cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells);
Timo Koch's avatar
Timo Koch committed
579

580
581
            // periodic boundaries
            const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>());
582

583
            // get the overlap dependent on the discretization method
584
            const int overlap = YaspOverlapHelper<discMethod>::getOverlap(modelParamGroup);
585

586
587
588
589
590
591
592
593
594
595
596
597
598
599
            // make the grid
            if (!haveParamInGroup(modelParamGroup, "Grid.Partitioning"))
            {
                // construct using default load balancing
                ParentType::gridPtr() = std::make_shared<Grid>(upperRight, cells, periodic, overlap);
            }
            else
            {
                // construct using user defined partitioning
                const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
                Dune::YaspFixedSizePartitioner<dim> lb(partitioning);
                ParentType::gridPtr() = std::make_shared<Grid>(upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
            }
            postProcessing_(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
600
        }
601
602

        // Didn't find a way to construct the grid
603
604
        else
        {
605
606
607
608
609
            const auto prefix = modelParamGroup == "" ? modelParamGroup : modelParamGroup + ".";
            DUNE_THROW(ParameterException, "Please supply one of the parameters "
                                           << prefix + "Grid.UpperRight"
                                           << ", or a grid file in " << prefix + "Grid.File");

Timo Koch's avatar
Timo Koch committed
610
611
612
613
614
615
616
        }
    }

private:
    /*!
     * \brief Postprocessing for YaspGrid
     */
617
    static void postProcessing_(const std::string& modelParamGroup)
Timo Koch's avatar
Timo Koch committed
618
619
    {
        // Check if should refine the grid
620
        bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true);
Timo Koch's avatar
Timo Koch committed
621
        ParentType::grid().refineOptions(keepPhysicalOverlap);
622
        ParentType::maybeRefineGrid(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
    }
};

/*!
 * \brief Provides a grid creator for YaspGrids with non-zero offset
 *        from information in the input file
 *
 * All keys are expected to be in group GridParameterGroup.
 * The following keys are recognized:
 * - LowerLeft : lower left corner coordinates
 * - UpperRight : upper right corner coordinates
 * - Cells : the number of cells in each direction
 * - Periodic : true or false for each direction
 * - Overlap : overlap size in cells
 * - Partitioning : a non-standard load-balancing, number of processors per direction
 * - KeepPyhsicalOverlap : whether to keep the physical overlap
 *     in physical size or in number of cells upon refinement
 * - Refinement : the number of global refines to apply initially.
 *
 */
643
644
template<DiscretizationMethod discMethod, class ct, int dim>
class GridCreatorImpl<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim> >, discMethod>
645
          : public GridCreatorBase<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim> > >
Timo Koch's avatar
Timo Koch committed
646
647
{
public:
648
649
    using Grid = typename Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim> >;
    using ParentType = GridCreatorBase<Grid>;
Timo Koch's avatar
Timo Koch committed
650
651
652
653

    /*!
     * \brief Make the grid. This is implemented by specializations of this method.
     */
654
    static void makeGrid(const std::string& modelParamGroup = "")
Timo Koch's avatar
Timo Koch committed
655
    {
656
657
658
659
660
661
662
663
664
        // First try to create it from a DGF file in GridParameterGroup.File
        if (haveParamInGroup(modelParamGroup, "Grid.File"))
        {
            ParentType::makeGridFromDgfFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"));
            postProcessing_(modelParamGroup);
            return;
        }

        // Then look for the necessary keys to construct from the input file
665
666
667
668
669
        else if (haveParamInGroup(modelParamGroup, "Grid.UpperRight"))
        {
            using GlobalPosition = Dune::FieldVector<ct, dim>;
            const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight");
            const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0));
670

671
672
673
            // number of cells in each direction
            std::array<int, dim> cells; cells.fill(1);
            cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells);
674

675
676
            // periodic boundaries
            const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>());
677

678
            // get the overlap dependent on some template parameters
679
            const int overlap = YaspOverlapHelper<discMethod>::getOverlap(modelParamGroup);
680

681
682
683
684
685
686
687
688
689
690
691
692
693
694
            // make the grid
            if (!haveParamInGroup(modelParamGroup, "Grid.Partitioning"))
            {
                // construct using default load balancing
                ParentType::gridPtr() = std::make_shared<Grid>(lowerLeft, upperRight, cells, periodic, overlap);
            }
            else
            {
                // construct using user defined partitioning
                const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
                Dune::YaspFixedSizePartitioner<dim> lb(partitioning);
                ParentType::gridPtr() = std::make_shared<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
            }
            postProcessing_(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
695
        }
696
697

        // Didn't find a way to construct the grid
698
699
        else
        {
700
701
702
703
704
            const auto prefix = modelParamGroup == "" ? modelParamGroup : modelParamGroup + ".";
            DUNE_THROW(ParameterException, "Please supply one of the parameters "
                                           << prefix + "Grid.UpperRight"
                                           << ", or a grid file in " << prefix + "Grid.File");

Timo Koch's avatar
Timo Koch committed
705
706
707
708
709
710
711
        }
    }

private:
    /*!
     * \brief Postprocessing for YaspGrid
     */
712
    static void postProcessing_(const std::string& modelParamGroup)
Timo Koch's avatar
Timo Koch committed
713
714
    {
        // Check if should refine the grid
715
        const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true);
Timo Koch's avatar
Timo Koch committed
716
        ParentType::grid().refineOptions(keepPhysicalOverlap);
717
        ParentType::maybeRefineGrid(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
718
719
720
    }
};

721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
/*!
 * \brief Provides a grid creator for YaspGrids with different zones and grading
 *
 * All keys are expected to be in group GridParameterGroup.
 * The following keys are recognized:
 * - Positions0 : position array for x-coordinate
 * - Positions1 : position array for y-coordinate
 * - Positions2 : position array for z-coordinate
 * - Cells0 : number of cells array for x-coordinate
 * - Cells1 : number of cells array for y-coordinate
 * - Cells2 : number of cells array for z-coordinate
 * - Grading0 : grading factor array for x-coordinate
 * - Grading1 : grading factor array for y-coordinate
 * - Grading2 : grading factor array for z-coordinate
 * - Verbosity : whether the grid construction should output to standard out
 * - Periodic : true or false for each direction
 * - Overlap : overlap size in cells
 * - Partitioning : a non-standard load-balancing, number of processors per direction
 * - KeepPyhsicalOverlap : whether to keep the physical overlap
 *     in physical size or in number of cells upon refinement
 * - Refinement : the number of global refines to apply initially.
 *
 * The grading factor \f$ g \f$ specifies the ratio between the next and the current cell size:
 * \f$ g = \frac{h_{i+1}}{h_i} \f$.
 * Negative grading factors are converted to
 * \f$ g = -\frac{1}{g_\textrm{negative}} \f$
 * to avoid issues with imprecise fraction numbers.
 */
749
750
template<DiscretizationMethod discMethod, class ctype, int dim>
class GridCreatorImpl<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >, discMethod>
751
          : public GridCreatorBase<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > >
752
753
{
public:
754
755
    using Grid = typename Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >;
    using ParentType = GridCreatorBase<Grid>;
756
757
758
759

    /*!
     * \brief Make the grid. This is implemented by specializations of this method.
     */
760
    static void makeGrid(const std::string& modelParamGroup = "")
761
762
763
    {
        // Only construction from the input file is possible
        // Look for the necessary keys to construct from the input file
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
        // The positions
        std::array<std::vector<ctype>, dim> positions;
        for (int i = 0; i < dim; ++i)
            positions[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Positions" + std::to_string(i));

        // the number of cells (has a default)
        std::array<std::vector<int>, dim> cells;
        for (int i = 0; i < dim; ++i)
        {
            cells[i].resize(positions[i].size()-1, 1.0);
            cells[i] = getParamFromGroup<std::vector<int>>(modelParamGroup, "Grid.Cells" + std::to_string(i), cells[i]);
        }

        // grading factor (has a default)
        std::array<std::vector<ctype>, dim> grading;
        for (int i = 0; i < dim; ++i)
        {
            grading[i].resize(positions[i].size()-1, 1.0);
            grading[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]);
        }

785
786
787
788
789
790
791
792
793
794
795
796
797
798
        // call the generic function
        makeGrid(positions, cells, grading, modelParamGroup);
    }

    /*!
     * \brief Make the grid using input data not read from the input file.
     */
    static void makeGrid(const std::array<std::vector<ctype>, dim>& positions,
                         const std::array<std::vector<int>, dim>& cells,
                         const std::array<std::vector<ctype>, dim>& grading,
                         const std::string& modelParamGroup = "")
    {


799
800
        // Additional arameters (they have a default)
        const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>());
801
        const int overlap = YaspOverlapHelper<discMethod>::getOverlap(modelParamGroup);
802
803
804
805
806
807
        const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);

        // Some sanity checks
        for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
        {
            if (cells[dimIdx].size() + 1 != positions[dimIdx].size())
808
            {
809
                DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Cells\" and \"Positions\" arrays");
810
            }
811
            if (grading[dimIdx].size() + 1 != positions[dimIdx].size())
812
            {
813
                DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Grading\" and \"Positions\" arrays");
814
            }
815
816
            ctype temp = std::numeric_limits<ctype>::lowest();
            for (unsigned int posIdx = 0; posIdx < positions[dimIdx].size(); ++posIdx)
817
            {
818
819
820
821
822
                if (temp > positions[dimIdx][posIdx])
                {
                    DUNE_THROW(Dune::RangeError, "Make sure to specify a monotone increasing \"Positions\" array");
                }
                temp = positions[dimIdx][posIdx];
823
            }
824
        }
825

826
        const auto globalPositions = computeGlobalPositions_(positions, cells, grading, verbose);
827

828
829
830
831
832
833
834
835
836
837
838
839
840
        // make the grid
        if (!haveParamInGroup(modelParamGroup, "Grid.Partitioning"))
        {
            // construct using default load balancing
            ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap);
        }
        else
        {
            // construct using user defined partitioning
            const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning");
            Dune::YaspFixedSizePartitioner<dim> lb(partitioning);
            ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap, typename Grid::CollectiveCommunicationType(), &lb);
        }
841

842
843
        postProcessing_(modelParamGroup);
    }
844

845
846
847
848
849
850
851
852
853
854
855
private:
    /*!
     * \brief Postprocessing for YaspGrid
     */
    static void postProcessing_(const std::string& modelParamGroup)
    {
        // Check if should refine the grid
        const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true);
        ParentType::grid().refineOptions(keepPhysicalOverlap);
        ParentType::maybeRefineGrid(modelParamGroup);
    }
856

857
858
859
860
861
862
863
864
865
866
867
868
    //! Compute the global position tensor grid from the given positions, cells, and grading factors
    static std::array<std::vector<ctype>, dim>
    computeGlobalPositions_(const std::array<std::vector<ctype>, dim>& positions,
                            const std::array<std::vector<int>, dim>& cells,
                            const std::array<std::vector<ctype>, dim>& grading,
                            bool verbose = false)
    {
        std::array<std::vector<ctype>, dim> globalPositions;
        using std::pow;
        for (int dimIdx = 0; dimIdx < dim; dimIdx++)
        {
            for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
869
            {
870
871
872
873
874
875
876
877
878
                ctype lower = positions[dimIdx][zoneIdx];
                ctype upper = positions[dimIdx][zoneIdx+1];
                int numCells = cells[dimIdx][zoneIdx];
                ctype gradingFactor = grading[dimIdx][zoneIdx];
                ctype length = upper - lower;
                ctype height = 1.0;
                bool increasingCellSize = false;

                if (verbose)
879
                {
880
881
882
883
884
                    std::cout << "dim " << dimIdx
                              << " lower "  << lower
                              << " upper "  << upper
                              << " numCells "  << numCells
                              << " grading "  << gradingFactor;
885
                }
886
887

                if (gradingFactor > 1.0)
888
                {
889
                    increasingCellSize = true;
890
                }
891
892
893
894

                // take absolute values and reverse cell size increment to achieve
                // reverse behavior for negative values
                if (gradingFactor < 0.0)
895
                {
896
897
898
                    using std::abs;
                    gradingFactor = abs(gradingFactor);
                    if (gradingFactor < 1.0)
899
                    {
900
                        increasingCellSize = true;
901
902
903
                    }
                }

904
905
                // if the grading factor is exactly 1.0 do equal spacing
                if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7)
906
                {
907
                    height = 1.0 / numCells;
908
909
                    if (verbose)
                    {
910
                        std::cout << " -> h "  << height * length << std::endl;
911
                    }
912
913
914
915
916
                }
                // if grading factor is not 1.0, do power law spacing
                else
                {
                    height = (1.0 - gradingFactor) / (1.0 - pow(gradingFactor, numCells));
917

918
                    if (verbose)
919
                    {
920
921
922
923
                        std::cout << " -> grading_eff "  << gradingFactor
                                  << " h_min "  << height * pow(gradingFactor, 0) * length
                                  << " h_max "  << height * pow(gradingFactor, numCells-1) * length
                                  << std::endl;
924
                    }
925
                }
926

927
928
929
930
931
932
                std::vector<ctype> localPositions;
                localPositions.push_back(0);
                for (int i = 0; i < numCells-1; i++)
                {
                    ctype hI = height;
                    if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7))
933
                    {
934
                        if (increasingCellSize)
935
                        {
936
                            hI *= pow(gradingFactor, i);
937
                        }
938
                        else
939
                        {
940
                            hI *= pow(gradingFactor, numCells-i-1);
941
942
                        }
                    }
943
944
                    localPositions.push_back(localPositions[i] + hI);
                }
945

946
947
948
949
                for (int i = 0; i < localPositions.size(); i++)
                {
                    localPositions[i] *= length;
                    localPositions[i] += lower;
950
951
                }

952
953
954
955
                for (unsigned int i = 0; i < localPositions.size(); ++i)
                {
                    globalPositions[dimIdx].push_back(localPositions[i]);
                }
956
            }
957
            globalPositions[dimIdx].push_back(positions[dimIdx].back());
958
959
        }

960
        return globalPositions;
961
962
963
    }
};

Timo Koch's avatar
Timo Koch committed
964
965
966
967
968
969
970
971
972
973
974
975
976
/*!
 * \brief Provides a grid creator for OneDGrids
 *        from information in the input file
 *
 * All keys are expected to be in group GridParameterGroup.
 * The following keys are recognized:
 * - LeftBoundary : start coordinate
 * - RightBoundary : end coordinate
 * - Cells : the number of cell
 * - RefinementType : local or copy
 * - Refinement : the number of global refines to apply initially.
 *
 */
977
978
template<DiscretizationMethod discMethod>
class GridCreatorImpl<Dune::OneDGrid, discMethod>
979
          : public GridCreatorBase<Dune::OneDGrid>
Timo Koch's avatar
Timo Koch committed
980
981
{
public:
982
983
    using Grid = Dune::OneDGrid;
    using ParentType = GridCreatorBase<Grid>;
Timo Koch's avatar
Timo Koch committed
984
985
986
987

    /*!
     * \brief Make the grid. This is implemented by specializations of this method.
     */
988
    static void makeGrid(const std::string& modelParamGroup = "")
Timo Koch's avatar
Timo Koch committed
989
    {
990
991
992
993
994
995

        // try to create it from a DGF or msh file in GridParameterGroup.File
        if (haveParamInGroup(modelParamGroup, "Grid.File"))
        {
            ParentType::makeGridFromDgfFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"));
            postProcessing_(modelParamGroup);
Timo Koch's avatar
Timo Koch committed
996
997
998
999
            return;
        }

        // Look for the necessary keys to construct from the input file
1000
        else if (haveParamInGroup(modelParamGroup, "Grid.RightBoundary"))
For faster browsing, not all history is shown. View entire blame