gridmanager_alu.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_IO_GRID_MANAGER_ALU_HH
#define DUMUX_IO_GRID_MANAGER_ALU_HH
 
// ALUGrid specific includes
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/dgf.hh>
#endif
 
#ifndef DUMUX_IO_GRID_MANAGER_BASE_HH
#include <dumux/io/grid/gridmanager_base.hh>
#endif
 
#include <dumux/common/boundaryflag.hh>
#include <dumux/common/gridcapabilities.hh>
 
namespace Dumux {
 
#if HAVE_DUNE_ALUGRID

template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refinementType>
class GridManager<Dune::ALUGrid<dim, dimworld, elType, refinementType>>
: public GridManagerBase<Dune::ALUGrid<dim, dimworld, elType, refinementType>>
{
public:
    using Grid = Dune::ALUGrid<dim, dimworld, elType, refinementType>;
    using ParentType = GridManagerBase<Grid>;

    void init(const std::string& modelParamGroup = "", bool adaptiveRestart = false)
    {
        // restarting an adaptive grid using Dune's BackupRestoreFacility
        // TODO: the part after first || is backward compatibility with old sequential models remove once sequential adaptive restart is replaced
        if (adaptiveRestart || hasParam("Restart") || hasParam("TimeManager.Restart"))
        {
            auto restartTime = getParamFromGroup<double>(modelParamGroup, "TimeLoop.Restart", 0.0);
            // TODO: backward compatibility with old sequential models remove once sequential adaptive restart is replaced
            if (hasParam("Restart") || hasParam("TimeManager.Restart"))
            {
                restartTime = getParamFromGroup<double>("TimeManager", "Restart");
                std::cerr << "Warning: You are using a deprecated restart mechanism. The usage will change in the future.\n";
            }
 
            const int rank = Dune::MPIHelper::getCommunication().rank();
            const std::string name = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name");
            std::ostringstream oss;
            oss << name << "_time=" << restartTime << "_rank=" << rank << ".grs";
            std::cout << "Restoring an ALUGrid from " << oss.str() << std::endl;
            ParentType::gridPtr() = std::shared_ptr<Grid>(Dune::BackupRestoreFacility<Grid>::restore(oss.str()));
            ParentType::loadBalance();
            return;
        }
 
        // First, try to create it from a file in GridParameterGroup.File
        else if (hasParamInGroup(modelParamGroup, "Grid.File"))
        {
            makeGridFromFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"), modelParamGroup);
            ParentType::maybeRefineGrid(modelParamGroup);
            ParentType::loadBalance();
            return;
        }
 
        // Then look for the necessary keys to construct from the input file
        else if (hasParamInGroup(modelParamGroup, "Grid.UpperRight"))
        {
            if (elType == Dune::cube)
                makeStructuredGrid<dim, dimworld>(ParentType::CellType::Cube, modelParamGroup);
            else if (elType == Dune::simplex)
                makeStructuredGrid<dim, dimworld>(ParentType::CellType::Simplex, modelParamGroup);
            else
                DUNE_THROW(Dune::IOError, "ALUGrid only supports Dune::cube or Dune::simplex as cell type!");
 
            ParentType::maybeRefineGrid(modelParamGroup);
            ParentType::loadBalance();
        }
 
        // Didn't find a way to construct the grid
        else
        {
            const auto prefix = modelParamGroup.empty() ? modelParamGroup : modelParamGroup + ".";
            DUNE_THROW(ParameterException, "Please supply one of the parameters "
                                           << prefix + "Grid.UpperRight"
                                           << ", or a grid file in " << prefix + "Grid.File");
 
        }
    }

    void makeGridFromFile(const std::string& fileName,
                          const std::string& modelParamGroup)
    {
        // We found a file in the input file...does it have a supported extension?
        const std::string extension = ParentType::getFileExtension(fileName);
        if (extension != "dgf" && extension != "msh" && extension != "vtu" && extension != "vti")
            DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " doesn't support grid files with extension: *."<< extension);
 
        // Dune Grid Format (DGF) files
        if (extension == "dgf")
        {
            ParentType::enableDgfGridPointer_ = true;
            ParentType::dgfGridPtr() = Dune::GridPtr<Grid>(fileName.c_str(), Dune::MPIHelper::getCommunicator());
            ParentType::gridData_ = std::make_shared<typename ParentType::GridData>(ParentType::dgfGridPtr());
        }
 
        // Gmsh mesh format
        else if (extension == "msh")
        {
            // get some optional parameters
            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);
 
            if (domainMarkers)
            {
                ParentType::enableGmshDomainMarkers_ = true;
                std::vector<int> boundaryMarkersInsertionIndex, boundaryMarkers, faceMarkers, elementMarkers;
                auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>();
                Dune::GmshReader<Grid>::read(*gridFactory, fileName, boundaryMarkersInsertionIndex, elementMarkers, verbose, boundarySegments);
                ParentType::gridPtr() = std::shared_ptr<Grid>(gridFactory->createGrid());
 
                // reorder boundary markers according to boundarySegmentIndex
                boundaryMarkers.resize(ParentType::gridPtr()->numBoundarySegments(), 0);
                faceMarkers.resize(ParentType::gridPtr()->leafGridView().size(1), 0);
                const auto& indexSet = ParentType::gridPtr()->leafGridView().indexSet();
                for (const auto& element : elements(ParentType::gridPtr()->leafGridView()))
                {
                    for (const auto& intersection : intersections(ParentType::gridPtr()->leafGridView(), element))
                    {
                        if (intersection.boundary() && gridFactory->wasInserted(intersection))
                        {
                            auto marker = boundaryMarkersInsertionIndex[gridFactory->insertionIndex(intersection)];
                            boundaryMarkers[intersection.boundarySegmentIndex()] = marker;
                            faceMarkers[indexSet.index(element.template subEntity<1>(intersection.indexInInside()))] = marker;
                        }
                    }
                }
 
                ParentType::gridData_ = std::make_shared<typename ParentType::GridData>(ParentType::gridPtr(), std::move(gridFactory),
                                                       std::move(elementMarkers), std::move(boundaryMarkers), std::move(faceMarkers));
            }
            else
            {
                auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>();
                Dune::GmshReader<Grid>::read(*gridFactory, fileName, verbose, boundarySegments);
                ParentType::gridPtr() = std::shared_ptr<Grid>(gridFactory->createGrid());
            }
        }
 
        // VTK file formats for unstructured grids
        // (can be constructed from both unstructured or structured grid data file formats)
        else if (extension == "vtu" || extension == "vti")
        {
            VTKReader vtkReader(fileName);
            VTKReader::Data cellData, pointData;
            auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>();
            const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
            ParentType::gridPtr() = vtkReader.readGrid(*gridFactory, cellData, pointData, verbose);
            ParentType::gridData_ = std::make_shared<typename ParentType::GridData>(ParentType::gridPtr(), std::move(gridFactory), std::move(cellData), std::move(pointData));
            ParentType::enableVtkData_ = true;
        }
    }

    template <int dimension, int dimensionworld, std::enable_if_t<dimension != dimensionworld, int> = 0>
    void makeStructuredGrid(typename ParentType::CellType cellType,
                            const std::string& modelParamGroup)
    {
        DUNE_THROW(Dune::IOError, "ALUGrid currently only supports the creation of structured grids with dimension == dimensionworld. Consider reading in a grid file instead.");
    }

    template <int dimension, int dimensionworld, std::enable_if_t<dimension == dimensionworld, int> = 0>
    void makeStructuredGrid(typename ParentType::CellType cellType,
                            const std::string& modelParamGroup)
    {
        // make a structured grid
        if (elType == Dune::cube)
            ParentType::template makeStructuredGrid<dimension, dimensionworld>(ParentType::CellType::Cube, modelParamGroup);
        else if (elType == Dune::simplex)
            ParentType::template makeStructuredGrid<dimension, dimensionworld>(ParentType::CellType::Simplex, modelParamGroup);
        else
            DUNE_THROW(Dune::IOError, "ALUGrid only supports Dune::cube or Dune::simplex as cell type!");
    }
};

template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refinementType>
class BoundaryFlag<Dune::ALUGrid<dim, dimworld, elType, refinementType>>
{
public:
    BoundaryFlag() : flag_(invalidFlag_) {}
 
    template<class Intersection>
    BoundaryFlag(const Intersection& i) : flag_(invalidFlag_)
    {
        if (i.boundary())
            flag_ = i.impl().boundaryId();
    }
 
    using value_type = int;
 
    value_type get() const { return flag_; }
 
    operator bool() const { return flag_ != invalidFlag_; }
 
private:
    static constexpr value_type invalidFlag_ = -1;
    value_type flag_;
};
 
namespace Grid::Capabilities {
 
// To the best of our knowledge ALUGrid is view thread-safe
// This specialization can be removed after we depend on Dune release 2.9 in which this is guaranteed by ALUGrid itself
template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refinementType>
struct MultithreadingSupported<Dune::ALUGrid<dim, dimworld, elType, refinementType>>
{
    template<class GV>
    static bool eval(const GV&) // default is independent of the grid view
    { return true; }
};
 
} // end namespace Grid::Capabilities
 
#endif // HAVE_DUNE_ALUGRID
 
} // end namespace Dumux
 
#endif