diff --git a/dumux/io/grid/gmshgriddatahandle.hh b/dumux/io/grid/gmshgriddatahandle.hh new file mode 100644 index 0000000000000000000000000000000000000000..a5a932e884c8815f1a92151be3400c165ec9c0f6 --- /dev/null +++ b/dumux/io/grid/gmshgriddatahandle.hh @@ -0,0 +1,188 @@ +// -*- 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 + * \ingroup InputOutput + * \brief A data handle for commucating grid data for gmsh grids + */ +#ifndef DUMUX_GMSH_GRID_DATA_HANDLE_HH +#define DUMUX_GMSH_GRID_DATA_HANDLE_HH + +#include <memory> +#include <algorithm> + +#include <dune/common/parallel/collectivecommunication.hh> +#include <dune/grid/common/partitionset.hh> +#include <dune/grid/common/datahandleif.hh> + +// UGGrid specific includes +#if HAVE_UG +#include <dune/grid/uggrid.hh> +#endif + +namespace Dumux { + +/*! + * \ingroup InputOutput + * \brief A data handle for commucating grid data for gmsh grids + */ +template<class Grid, class GridFactory, class Data> +struct GmshGridDataHandle : public Dune::CommDataHandleIF<GmshGridDataHandle<Grid, GridFactory, Data>, typename Data::value_type> +{ + using GridView = typename Grid::LevelGridView; + + GmshGridDataHandle(const Grid& grid, const GridFactory& gridFactory, Data& elementMarkers, Data& boundaryMarkers, Data& faceMarkers) + : gridView_(grid.levelGridView(0)) + , idSet_(grid.localIdSet()) + , elementMarkers_(elementMarkers) + , boundaryMarkers_(boundaryMarkers) + , faceMarkers_(faceMarkers) + { + const auto& indexSet = gridView_.indexSet(); + + for (const auto& element : elements(gridView_, Dune::Partitions::interior)) + std::swap(elementMarkers_[gridFactory.insertionIndex(element)], data_[idSet_.id(element)]); + + for (const auto& face : entities(gridView_, Dune::Codim<1>())) + std::swap(faceMarkers_[indexSet.index(face)], data_[idSet_.id(face)]); + } + + ~GmshGridDataHandle() + { + const auto& indexSet = gridView_.indexSet(); + + elementMarkers_.resize(indexSet.size(0)); + for (const auto& element : elements(gridView_)) + std::swap(elementMarkers_[indexSet.index(element)], data_[idSet_.id(element)]); + + faceMarkers_.resize(indexSet.size(1)); + for (const auto& face : entities(gridView_, Dune::Codim<1>())) + std::swap(faceMarkers_[indexSet.index(face)], data_[idSet_.id(face)]); + + boundaryMarkers_.resize(gridView_.grid().numBoundarySegments(), 0); + for (const auto& element : elements(gridView_.grid().leafGridView())) + { + for (const auto& intersection : intersections(gridView_.grid().leafGridView(), element)) + { + if (intersection.boundary()) + { + const auto marker = faceMarkers_[indexSet.index(element.template subEntity<1>(intersection.indexInInside()))]; + boundaryMarkers_[intersection.boundarySegmentIndex()] = marker; + } + } + } + } + + Dune::CommDataHandleIF<GmshGridDataHandle<Grid, GridFactory, Data>, typename Data::value_type>& interface() + { return *this; } + + bool contains (int dim, int codim) const + { return codim == 0 || codim == 1; } + + bool fixedSize (int dim, int codim) const + { return true; } + + template<class EntityType> + std::size_t size (const EntityType& e) const + { return 1; } + + template<class MessageBufferImp, class EntityType> + void gather (MessageBufferImp& buff, const EntityType& e) const + { buff.write(data_[idSet_.id(e)]); } + + template<class MessageBufferImp, class EntityType> + void scatter (MessageBufferImp& buff, const EntityType& e, std::size_t n) + { buff.read(data_[idSet_.id(e)]); } + +private: + using IdSet = typename Grid::LocalIdSet; + + const GridView gridView_; + const IdSet &idSet_; + Data& elementMarkers_; + Data& boundaryMarkers_; + Data& faceMarkers_; + mutable std::map< typename IdSet::IdType, typename Data::value_type> data_; +}; + +#if HAVE_UG + +/*! + * \ingroup InputOutput + * \brief A data handle for commucating grid data for gmsh grids (specialization for UGGrid) + */ +template<class GridFactory, class Data, int dimgrid> +struct GmshGridDataHandle<Dune::UGGrid<dimgrid>, GridFactory, Data> +: public Dune::CommDataHandleIF<GmshGridDataHandle<Dune::UGGrid<dimgrid>, GridFactory, Data>, typename Data::value_type> +{ + using Grid = Dune::UGGrid<dimgrid>; + using GridView = typename Grid::LevelGridView; + + GmshGridDataHandle(const Grid& grid, const GridFactory& gridFactory, Data& elementMarkers) + : gridView_(grid.levelGridView(0)) + , idSet_(grid.localIdSet()) + , elementMarkers_(elementMarkers) + { + for (const auto& element : elements(gridView_, Dune::Partitions::interior)) + std::swap(elementMarkers_[gridFactory.insertionIndex(element)], data_[idSet_.id(element)]); + } + + ~GmshGridDataHandle() + { + const auto& indexSet = gridView_.indexSet(); + elementMarkers_.resize(indexSet.size(0)); + for (const auto& element : elements(gridView_)) + std::swap(elementMarkers_[indexSet.index(element)], data_[idSet_.id(element)]); + } + + Dune::CommDataHandleIF<GmshGridDataHandle<Grid, GridFactory, Data>, typename Data::value_type>& interface() + { return *this; } + + bool contains (int dim, int codim) const + { return codim == 0 || codim == 1; } + + bool fixedSize (int dim, int codim) const + { return true; } + + template<class EntityType> + std::size_t size (const EntityType& e) const + { return 1; } + + template<class MessageBufferImp, class EntityType> + void gather (MessageBufferImp& buff, const EntityType& e) const + { buff.write(data_[idSet_.id(e)]); } + + template<class MessageBufferImp, class EntityType> + void scatter (MessageBufferImp& buff, const EntityType& e, std::size_t n) + { buff.read(data_[idSet_.id(e)]); } + +private: + using IdSet = typename Grid::LocalIdSet; + + const GridView gridView_; + const IdSet &idSet_; + Data& elementMarkers_; + mutable std::map< typename IdSet::IdType, typename Data::value_type> data_; +}; + +#endif // HAVE_UG + +} // namespace Dumux + +#endif diff --git a/dumux/io/grid/griddata.hh b/dumux/io/grid/griddata.hh new file mode 100644 index 0000000000000000000000000000000000000000..d7a7556f631e19941e81778ba25af3b189f14353 --- /dev/null +++ b/dumux/io/grid/griddata.hh @@ -0,0 +1,210 @@ +// -*- 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 + * \ingroup InputOutput + * \brief Class for grid data attached to dgf or gmsh grid files + */ +#ifndef DUMUX_IO_GRID_DATA_HH +#define DUMUX_IO_GRID_DATA_HH + +#include <vector> +#include <memory> +#include <type_traits> + +#include <dune/common/exceptions.hh> +#include <dune/common/parallel/collectivecommunication.hh> +#include <dune/common/parallel/mpihelper.hh> +#include <dune/grid/common/gridfactory.hh> + +// UGGrid specific includes +#if HAVE_UG +#include <dune/grid/uggrid.hh> +#endif + +#include "gmshgriddatahandle.hh" + +namespace Dumux { + +namespace Detail { + +template<class Grid> +struct isUG : public std::false_type {}; + +#if HAVE_UG +template<int dim> +struct isUG<Dune::UGGrid<dim>> : public std::true_type {}; +#endif + +} // end namespace Details + +/*! + * \ingroup InputOutput + * \brief Class for grid data attached to dgf or gmsh grid files + */ +template <class Grid> +class GridData +{ + using Intersection = typename Grid::LeafIntersection; + using Element = typename Grid::template Codim<0>::Entity; + using DataHandle = GmshGridDataHandle<Grid, Dune::GridFactory<Grid>, std::vector<int>>; + +public: + //! constructor for gmsh grid data + GridData(std::shared_ptr<Grid> grid, std::shared_ptr<Dune::GridFactory<Grid>> factory, + std::vector<int>&& elementMarkers, std::vector<int>&& boundaryMarkers, std::vector<int>&& faceMarkers = std::vector<int>{}) + : gmshGrid_(grid) + , gridFactory_(factory) + , elementMarkers_(elementMarkers) + , boundaryMarkers_(boundaryMarkers) + , faceMarkers_(faceMarkers) + {} + + //! constructor for dgf grid data + GridData(Dune::GridPtr<Grid> grid) + : dgfGrid_(grid) + , isDgfData_(true) + {} + + + /*! + * \brief Call the parameters function of the DGF grid pointer if available + */ + template <class Entity> + const std::vector<double>& parameters(const Entity& entity) const + { + if (isDgfData_) + { + if (entity.hasFather()) + { + auto level0entity = entity; + while(level0entity.hasFather()) + level0entity = level0entity.father(); + + + return dgfGrid_.parameters(level0entity); + } + else + { + return dgfGrid_.parameters(entity); + } + } + 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> + const Dune::DGFBoundaryParameter::type& parameters(const Dune::Intersection<GridImp, IntersectionImp>& intersection) const + { + if (isDgfData_) + return dgfGrid_.parameters(intersection); + else + DUNE_THROW(Dune::InvalidStateException, "The parameters method is only available if the grid was constructed with a DGF file."); + } + + /*! + * \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() + */ + int getBoundaryDomainMarker(int boundarySegmentIndex) const + { + if (!gmshGrid_) + DUNE_THROW(Dune::InvalidStateException, "Domain markers are only available for gmsh grids."); + if (boundarySegmentIndex >= boundaryMarkers_.size()) + DUNE_THROW(Dune::RangeError, "Boundary segment index "<< boundarySegmentIndex << " bigger than number of boundary segments in grid."); + return boundaryMarkers_[boundarySegmentIndex]; + } + + /*! + * \brief Return the boundary domain marker (Gmsh physical entity number) of an intersection + Only available when using Gmsh with GridParameterGroup.DomainMarkers = 1. + * \param intersection The intersection to be evaluated + */ + int getBoundaryDomainMarker(const Intersection& intersection) const + { return getBoundaryDomainMarker(intersection.boundarySegmentIndex()); } + + + /*! + * \brief Return the element domain marker (Gmsh physical entity number) of an element. + Only available when using Gmsh with GridParameterGroup.DomainMarkers = 1. + * \param element The element to be evaluated + */ + int getElementDomainMarker(const Element& element) const + { + if (!gmshGrid_) + DUNE_THROW(Dune::InvalidStateException, "Domain markers are only available for gmsh grids."); + + // parameters are only given for level 0 elements + 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 + // for UGGrid element data is read on all processes since UGGrid can't communicate element data (yet) + if (gmshGrid_->comm().size() > 1 && !Detail::isUG<Grid>::value) + return elementMarkers_[gmshGrid_->levelGridView(0).indexSet().index(level0element)]; + else + return elementMarkers_[gridFactory_->insertionIndex(level0element)]; + } + + /*! + * \brief Create a data handle for communication of the data in parallel simulations + * \note this data hande is the default + */ + template<bool ug = Detail::isUG<Grid>::value, typename std::enable_if_t<!ug, int> = 0> + DataHandle createGmshDataHandle() + { + return DataHandle(*gmshGrid_, *gridFactory_, elementMarkers_, boundaryMarkers_, faceMarkers_); + } + + /*! + * \brief Create a data handle for communication of the data in parallel simulations + * \note this data hande is the specialized for UGGrid since UGGrid can't communicate element data (yet) + */ + template<bool ug = Detail::isUG<Grid>::value, typename std::enable_if_t<ug, int> = 0> + DataHandle createGmshDataHandle() + { + return DataHandle(*gmshGrid_, *gridFactory_, elementMarkers_); + } + +private: + // grid and grid factor for gmsh grid data + std::shared_ptr<Grid> gmshGrid_; + std::shared_ptr<Dune::GridFactory<Grid>> gridFactory_; + + /*! + * \brief Element and boundary domain markers obtained from Gmsh physical entities + * They map from element indices / boundary ids to the physical entity number + */ + std::vector<int> elementMarkers_; + std::vector<int> boundaryMarkers_; + std::vector<int> faceMarkers_; + + // dgf grid data + Dune::GridPtr<Grid> dgfGrid_; + bool isDgfData_ = false; +}; + +} // namespace Dumux + +#endif diff --git a/dumux/io/grid/gridmanager.hh b/dumux/io/grid/gridmanager.hh new file mode 100644 index 0000000000000000000000000000000000000000..c9118be37edb40cc637ab0162680d86913ec5800 --- /dev/null +++ b/dumux/io/grid/gridmanager.hh @@ -0,0 +1,1285 @@ +// -*- 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 + * \ingroup InputOutput + * \brief Provides a grid manager for all supported grid managers with + * input file interfaces. Manages data via the grid data member. + * + * \todo add Petrel grids with dune-cornerpoint + */ +#ifndef DUMUX_IO_GRID_MANAGER_HH +#define DUMUX_IO_GRID_MANAGER_HH + +#include <array> +#include <bitset> +#include <memory> +#include <sstream> + +#include <dune/common/exceptions.hh> +#include <dune/common/classname.hh> +#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> +#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> +#include <dune/foamgrid/dgffoam.hh> +#endif + +#include <dumux/common/parameters.hh> +#include <dumux/discretization/methods.hh> + +#include "griddata.hh" + +namespace Dumux { + +/*! + * \ingroup InputOutput + * \brief The grid manager base interface (public) and methods common + * to most grid manager specializations (protected). + */ +template <class GridType> +class GridManagerBase +{ +public: + using Grid = GridType; + using GridData = Dumux::GridData<Grid>; + + /*! + * \brief Make the grid. Implement this method in the specialization of this class for a grid type. + */ + void init(const std::string& modelParamGroup = "") + { + DUNE_THROW(Dune::NotImplemented, + "The GridManager for grid type " << Dune::className<Grid>() << " is not implemented! Consider providing your own GridManager."); + } + + /*! + * \brief Returns a reference to the grid. + */ + Grid& grid() + { + if(enableDgfGridPointer_) + return *dgfGridPtr(); + else + return *gridPtr(); + } + + /*! + * \brief Call loadBalance() function of the grid. + */ + void loadBalance() + { + if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) + { + // if we may have dgf parameters use load balancing of the dgf pointer + if(enableDgfGridPointer_) + { + dgfGridPtr().loadBalance(); + // update the grid data + gridData_ = std::make_shared<GridData>(dgfGridPtr()); + } + + // if we have gmsh parameters we have to manually load balance the data + else if (enableGmshDomainMarkers_) + { + // element and face markers are communicated during load balance + auto dh = gridData_->createGmshDataHandle(); + gridPtr()->loadBalance(dh.interface()); + gridPtr()->communicate(dh.interface(), Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication); + } + else + gridPtr()->loadBalance(); + } + } + + std::shared_ptr<GridData> getGridData() const + { + if (!enableDgfGridPointer_ && !enableGmshDomainMarkers_) + DUNE_THROW(Dune::IOError, "No grid data available"); + + return gridData_; + } + + +protected: + + /*! + * \brief Returns a reference to the grid pointer (std::shared_ptr<Grid>) + */ + std::shared_ptr<Grid>& gridPtr() + { + if(!enableDgfGridPointer_) + return gridPtr_; + else + DUNE_THROW(Dune::InvalidStateException, "You are using DGF. To get the grid pointer use method dgfGridPtr()!"); + } + + /*! + * \brief Returns a reference to the DGF grid pointer (Dune::GridPtr<Grid>). + */ + Dune::GridPtr<Grid>& dgfGridPtr() + { + if(enableDgfGridPointer_) + return dgfGridPtr_; + else + DUNE_THROW(Dune::InvalidStateException, "The DGF grid pointer is only available if the grid was constructed with a DGF file!"); + } + + /*! + * \brief Returns the filename extension of a given filename + */ + std::string getFileExtension(const std::string& fileName) const + { + 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). + */ + 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 = getFileExtension(fileName); + if (extension != "dgf" && extension != "msh") + 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); + + // make the grid + if (extension == "dgf") + { + enableDgfGridPointer_ = true; + dgfGridPtr() = Dune::GridPtr<Grid>(fileName.c_str(), Dune::MPIHelper::getCommunicator()); + gridData_ = std::make_shared<GridData>(dgfGridPtr_); + } + 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) + enableGmshDomainMarkers_ = true; + + // as default read it on all processes in parallel + if(domainMarkers) + { + std::vector<int> boundaryMarkers, elementMarkers; + auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>(); + Dune::GmshReader<Grid>::read(*gridFactory, fileName, boundaryMarkers, elementMarkers, verbose, boundarySegments); + gridPtr() = std::shared_ptr<Grid>(gridFactory->createGrid()); + gridData_ = std::make_shared<GridData>(gridPtr_, std::move(gridFactory), std::move(elementMarkers), std::move(boundaryMarkers)); + } + else + { + auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>(); + Dune::GmshReader<Grid>::read(*gridFactory, fileName, verbose, boundarySegments); + gridPtr() = std::shared_ptr<Grid>(gridFactory->createGrid()); + } + } + } + + /*! + * \brief Makes a grid from a DGF file. This is used by grid managers that only support DGF. + */ + void makeGridFromDgfFile(const std::string& fileName) + { + // We found a file in the input file...does it have a supported extension? + const std::string extension = getFileExtension(fileName); + if(extension != "dgf") + DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " only supports DGF (*.dgf) but the specified filename has extension: *."<< extension); + + 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> + void makeStructuredGrid(CellType cellType, + const std::string& modelParamGroup) + { + 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)); + + using CellArray = std::array<unsigned int, dim>; + CellArray cells; cells.fill(1); + cells = getParamFromGroup<CellArray>(modelParamGroup, "Grid.Cells", cells); + + // make the grid + if (cellType == CellType::Cube) + { + gridPtr() = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, cells); + } + else if (cellType == CellType::Simplex) + { + gridPtr() = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, cells); + } + else + { + DUNE_THROW(Dune::GridError, "Unknown cell type for making structured grid! Choose Cube or Simplex."); + } + } + + /*! + * \brief Refines a grid after construction if GridParameterGroup.Refinement is set in the input file + */ + void maybeRefineGrid(const std::string& modelParamGroup) + { + if (haveParamInGroup(modelParamGroup, "Grid.Refinement")) + grid().globalRefine(getParamFromGroup<int>(modelParamGroup, "Grid.Refinement")); + } + + /*! + * \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. + */ + bool enableDgfGridPointer_ = false; + + /*! + * \brief A state variable if domain markers have been read from a Gmsh file. + */ + bool enableGmshDomainMarkers_ = false; + + std::shared_ptr<Grid> gridPtr_; + Dune::GridPtr<Grid> dgfGridPtr_; + + std::shared_ptr<GridData> gridData_; +}; + +/*! + * \brief The grid manager (this is the class used by the user) for all supported grid managers that constructs a grid + * from information in the input file and handles the data. + * \note This class is specialised below for all supported grid managers. It inherits the functionality of the base class. + */ +template <class Grid> +class GridManager +: public GridManagerBase<Grid> +{}; + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Specializations ////////////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/*! + * \brief Provides a grid manager 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. + * + */ +template<class ct, int dim> +class GridManager<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> >> +: public GridManagerBase<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> > > +{ +public: + using Grid = typename Dune::YaspGrid<dim, Dune::EquidistantCoordinates<ct, dim> >; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + // 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 + else if (haveParamInGroup(modelParamGroup, "Grid.UpperRight")) + { + + // get the upper right corner coordinates + const auto upperRight = getParamFromGroup<Dune::FieldVector<ct, dim>>(modelParamGroup, "Grid.UpperRight"); + + // number of cells in each direction + std::array<int, dim> cells; cells.fill(1); + cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells); + + // periodic boundaries + const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>()); + + // get the overlap + const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); + + // 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); + } + + // Didn't find a way to construct the grid + else + { + 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"); + + } + } + +private: + /*! + * \brief Postprocessing for YaspGrid + */ + void postProcessing_(const std::string& modelParamGroup) + { + // Check if should refine the grid + bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true); + ParentType::grid().refineOptions(keepPhysicalOverlap); + ParentType::maybeRefineGrid(modelParamGroup); + ParentType::loadBalance(); + } +}; + +/*! + * \brief Provides a grid manager 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. + * + */ +template<class ct, int dim> +class GridManager<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim>>> +: public GridManagerBase<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim>>> +{ +public: + using Grid = typename Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<ct, dim>>; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + // 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 + 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)); + + // number of cells in each direction + std::array<int, dim> cells; cells.fill(1); + cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells); + + // periodic boundaries + const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>()); + + // get the overlap dependent on some template parameters + const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); + + // 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); + } + + // Didn't find a way to construct the grid + else + { + 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"); + + } + } + +private: + /*! + * \brief Postprocessing for YaspGrid + */ + 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); + ParentType::loadBalance(); + } +}; + +/*! + * \brief Provides a grid manager 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. + */ +template<class ctype, int dim> +class GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >> +: public GridManagerBase<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > > +{ +public: + using Grid = typename Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + // Only construction from the input file is possible + // Look for the necessary keys to construct from the input file + // 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]); + } + + // call the generic function + init(positions, cells, grading, modelParamGroup); + } + + /*! + * \brief Make the grid using input data not read from the input file. + */ + void init(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 = "") + { + + + // Additional arameters (they have a default) + const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>()); + const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); + 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()) + { + DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Cells\" and \"Positions\" arrays"); + } + if (grading[dimIdx].size() + 1 != positions[dimIdx].size()) + { + DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Grading\" and \"Positions\" arrays"); + } + ctype temp = std::numeric_limits<ctype>::lowest(); + for (unsigned int posIdx = 0; posIdx < positions[dimIdx].size(); ++posIdx) + { + if (temp > positions[dimIdx][posIdx]) + { + DUNE_THROW(Dune::RangeError, "Make sure to specify a monotone increasing \"Positions\" array"); + } + temp = positions[dimIdx][posIdx]; + } + } + + const auto globalPositions = computeGlobalPositions_(positions, cells, grading, verbose); + + // 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); + } + + postProcessing_(modelParamGroup); + } + +private: + /*! + * \brief Postprocessing for YaspGrid + */ + 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); + ParentType::loadBalance(); + } + + //! Compute the global position tensor grid from the given positions, cells, and grading factors + 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) + { + 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) + { + std::cout << "dim " << dimIdx + << " lower " << lower + << " upper " << upper + << " numCells " << numCells + << " grading " << gradingFactor; + } + + if (gradingFactor > 1.0) + { + increasingCellSize = true; + } + + // take absolute values and reverse cell size increment to achieve + // reverse behavior for negative values + if (gradingFactor < 0.0) + { + using std::abs; + gradingFactor = abs(gradingFactor); + if (gradingFactor < 1.0) + { + increasingCellSize = true; + } + } + + // if the grading factor is exactly 1.0 do equal spacing + if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7) + { + height = 1.0 / numCells; + if (verbose) + { + std::cout << " -> h " << height * length << std::endl; + } + } + // if grading factor is not 1.0, do power law spacing + else + { + height = (1.0 - gradingFactor) / (1.0 - pow(gradingFactor, numCells)); + + if (verbose) + { + std::cout << " -> grading_eff " << gradingFactor + << " h_min " << height * pow(gradingFactor, 0) * length + << " h_max " << height * pow(gradingFactor, numCells-1) * length + << std::endl; + } + } + + 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)) + { + if (increasingCellSize) + { + hI *= pow(gradingFactor, i); + } + else + { + hI *= pow(gradingFactor, numCells-i-1); + } + } + localPositions.push_back(localPositions[i] + hI); + } + + for (int i = 0; i < localPositions.size(); i++) + { + localPositions[i] *= length; + localPositions[i] += lower; + } + + for (unsigned int i = 0; i < localPositions.size(); ++i) + { + globalPositions[dimIdx].push_back(localPositions[i]); + } + } + globalPositions[dimIdx].push_back(positions[dimIdx].back()); + } + + return globalPositions; + } +}; + +/*! + * \brief Provides a grid manager 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. + * + */ +template<> +class GridManager<Dune::OneDGrid> +: public GridManagerBase<Dune::OneDGrid> +{ +public: + using Grid = Dune::OneDGrid; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + + // 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); + return; + } + + // Look for the necessary keys to construct from the input file + else if (haveParamInGroup(modelParamGroup, "Grid.RightBoundary")) + { + // The required parameters + using CoordinateType = typename Grid::ctype; + const auto leftBoundary = getParamFromGroup<CoordinateType>(modelParamGroup, "Grid.LeftBoundary", 0.0); + const auto rightBoundary = getParamFromGroup<CoordinateType>(modelParamGroup, "Grid.RightBoundary"); + const int cells = getParamFromGroup<int>(modelParamGroup, "Grid.Cells", 1); + + ParentType::gridPtr() = std::make_shared<Grid>(cells, leftBoundary, rightBoundary); + postProcessing_(modelParamGroup); + return; + } + + // Look for the necessary keys to construct from the input file with just a coordinates vector + else if (haveParamInGroup(modelParamGroup, "Grid.Coordinates")) + { + const auto coordinates = getParamFromGroup<std::vector<typename Grid::ctype>>(modelParamGroup, "Grid.Coordinates"); + ParentType::gridPtr() = std::make_shared<Grid>(coordinates); + postProcessing_(modelParamGroup); + } + + // Didn't find a way to construct the grid + else + { + const auto prefix = modelParamGroup == "" ? modelParamGroup : modelParamGroup + "."; + DUNE_THROW(ParameterException, "Please supply one of the parameters " + << prefix + "Grid.RightBoundary" + << ", or " << prefix + "Grid.Coordinates" + << ", or a grid file in " << prefix + "Grid.File"); + } + } + + /*! + * \brief Call loadBalance() function of the grid. OneDGrid is not parallel an thus cannot communicate. + */ + void loadBalance() {} + +private: + /*! + * \brief Do some operatrion after making the grid, like global refinement + */ + void postProcessing_(const std::string& modelParamGroup) + { + // Set refinement type + const auto refType = getParamFromGroup<std::string>(modelParamGroup, "Grid.RefinementType", "Local"); + if (refType == "Local") + ParentType::grid().setRefinementType(Dune::OneDGrid::RefinementType::LOCAL); + else if (refType == "Copy") + ParentType::grid().setRefinementType(Dune::OneDGrid::RefinementType::COPY); + else + DUNE_THROW(Dune::IOError, "OneGrid only supports 'Local' or 'Copy' as refinment type. Not '"<< refType<<"'!"); + + // Check if should refine the grid + ParentType::maybeRefineGrid(modelParamGroup); + loadBalance(); + } +}; + +#if HAVE_UG + +/*! + * \brief Provides a grid manager for UGGrids + * from information in the input file + * + * All keys are expected to be in group GridParameterGroup. + + * The following keys are recognized: + * - File : A DGF or gmsh file to load from, type detection by file extension + * - LowerLeft : lowerLeft corner of a structured grid + * - UpperRight : upperright corner of a structured grid + * - Cells : number of elements in a structured grid + * - CellType : "Cube" or "Simplex" to be used for structured grids + * - Refinement : the number of global refines to perform + * - Verbosity : whether the grid construction should output to standard out + * - HeapSize: The heapsize used to allocate memory + * - BoundarySegments : whether to insert boundary segments into the grid + * + */ +template<int dim> +class GridManager<Dune::UGGrid<dim>> +: public GridManagerBase<Dune::UGGrid<dim>> +{ +public: + using Grid = typename Dune::UGGrid<dim>; + using ParentType = GridManagerBase<Grid>; + using Element = typename Grid::template Codim<0>::Entity; + + /*! + * \brief Make the UGGrid. + */ + void init(const std::string& modelParamGroup = "") + { + + // try to create it from a DGF or msh file in GridParameterGroup.File + if (haveParamInGroup(modelParamGroup, "Grid.File")) + { + preProcessing_(modelParamGroup); + ParentType::makeGridFromFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"), modelParamGroup); + postProcessing_(modelParamGroup); + return; + } + + // Then look for the necessary keys to construct from the input file + else if (haveParamInGroup(modelParamGroup, "Grid.UpperRight")) + { + preProcessing_(modelParamGroup); + // make the grid + const auto cellType = getParamFromGroup<std::string>(modelParamGroup, "Grid.CellType", "Cube"); + if (cellType == "Cube") + ParentType::template makeStructuredGrid<dim, dim>(ParentType::CellType::Cube, modelParamGroup); + else if (cellType == "Simplex") + ParentType::template makeStructuredGrid<dim, dim>(ParentType::CellType::Simplex, modelParamGroup); + else + DUNE_THROW(Dune::IOError, "UGGrid only supports 'Cube' or 'Simplex' as cell type. Not '"<< cellType<<"'!"); + postProcessing_(modelParamGroup); + } + + // Didn't find a way to construct the grid + else + { + 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"); + + } + } + + /*! + * \brief Call loadBalance() function of the grid. + * \note This overload is necessary because UGGrid cannot load balance element data yet. + * For dgf grids using element data in parallel will (hopefully) through an error + * For gmsh grids the parameters are read on every process so we use too much memory but + * the parameters are available via the insertionIndex of the level 0 element. + */ + void loadBalance() + { + if (Dune::MPIHelper::getCollectiveCommunication().size() > 1) + { + // if we may have dgf parameters use load balancing of the dgf pointer + if(ParentType::enableDgfGridPointer_) + { + ParentType::dgfGridPtr().loadBalance(); + // update the grid data + ParentType::gridData_ = std::make_shared<typename ParentType::GridData>(ParentType::dgfGridPtr()); + } + + // if we have gmsh parameters we have to manually load balance the data + else if (ParentType::enableGmshDomainMarkers_) + { + // element and face markers are communicated during load balance + auto dh = ParentType::gridData_->createGmshDataHandle(); + ParentType::gridPtr()->loadBalance(dh.interface()); + // Right now, UGGrid cannot communicate element data. If this gets implemented, communicate the data here: + // ParentType::gridPtr()->communicate(dh.interface(), Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication); + } + else + ParentType::gridPtr()->loadBalance(); + } + } + +private: + /*! + * \brief Do some operations before making the grid + */ + void preProcessing_(const std::string& modelParamGroup) + { + if(haveParamInGroup(modelParamGroup, "Grid.HeapSize")) + Grid::setDefaultHeapSize(getParamFromGroup<unsigned>(modelParamGroup, "Grid.HeapSize")); + } + + /*! + * \brief Do some operatrion after making the grid, like global refinement + */ + void postProcessing_(const std::string& modelParamGroup) + { + // Set refinement type + const auto refType = getParamFromGroup<std::string>(modelParamGroup, "Grid.RefinementType", "Local"); + if (refType == "Local") + ParentType::grid().setRefinementType(Dune::UGGrid<dim>::RefinementType::LOCAL); + else if (refType == "Copy") + ParentType::grid().setRefinementType(Dune::UGGrid<dim>::RefinementType::COPY); + else + DUNE_THROW(Dune::IOError, "UGGrid only supports 'Local' or 'Copy' as refinment type. Not '"<< refType<<"'!"); + + // Set closure type + const auto closureType = getParamFromGroup<std::string>(modelParamGroup, "Grid.ClosureType", "Green"); + if (closureType == "Green") + ParentType::grid().setClosureType(Dune::UGGrid<dim>::ClosureType::GREEN); + else if (closureType == "None") + ParentType::grid().setClosureType(Dune::UGGrid<dim>::ClosureType::NONE); + else + DUNE_THROW(Dune::IOError, "UGGrid only supports 'Green' or 'None' as closure type. Not '"<< closureType<<"'!"); + + // Check if should refine the grid + ParentType::maybeRefineGrid(modelParamGroup); + // do load balancing + loadBalance(); + } +}; + +#endif // HAVE_UG +// +#if HAVE_DUNE_ALUGRID + +/*! + * \brief Provides a grid manager for Dune ALUGrids + * from information in the input file + * + * All keys are expected to be in group GridParameterGroup. + + * The following keys are recognized: + * - File : A DGF or gmsh file to load from, type detection by file extension + * - LowerLeft : lowerLeft corner of a structured grid + * - UpperRight : upperright corner of a structured grid + * - Cells : number of elements in a structured grid + * - Refinement : the number of global refines to perform + * - Verbosity : whether the grid construction should output to standard out + * - BoundarySegments : whether to insert boundary segments into the grid + * + */ +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>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "", bool adaptiveRestart = false) + { + // restarting an adaptive grid using Dune's BackupRestoreFacility + // TODO: the part after first || is backward compatibilty with old sequential models remove once sequential adpative restart is replaced + if (adaptiveRestart || haveParam("Restart") || haveParam("TimeManager.Restart")) + { + auto restartTime = getParamFromGroup<double>(modelParamGroup, "TimeLoop.Restart", 0.0); + // TODO: backward compatibilty with old sequential models remove once sequential adpative restart is replaced + if (haveParam("Restart") || haveParam("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::getCollectiveCommunication().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; + } + + // try to create it from a DGF or msh file in GridParameterGroup.File + else if (haveParamInGroup(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 (haveParamInGroup(modelParamGroup, "Grid.UpperRight")) + { + // make a structured grid + if (elType == Dune::cube) + ParentType::template makeStructuredGrid<dim, dimworld>(ParentType::CellType::Cube, modelParamGroup); + else if (elType == Dune::simplex) + ParentType::template 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 == "" ? modelParamGroup : modelParamGroup + "."; + DUNE_THROW(ParameterException, "Please supply one of the parameters " + << prefix + "Grid.UpperRight" + << ", or a grid file in " << prefix + "Grid.File"); + + } + } + + /*! + * \brief Makes a grid from a file. We currently support *.dgf (Dune Grid Format) and *.msh (Gmsh mesh format). + */ + 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") + 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); + + // make the grid + 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()); + } + 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; + + // only filll the factory for rank 0 + if (domainMarkers) + { + std::vector<int> boundaryMarkersInsertionIndex, boundaryMarkers, faceMarkers, elementMarkers; + auto gridFactory = std::make_unique<Dune::GridFactory<Grid>>(); + if (Dune::MPIHelper::getCollectiveCommunication().rank() == 0) + 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>>(); + if (Dune::MPIHelper::getCollectiveCommunication().rank() == 0) + Dune::GmshReader<Grid>::read(*gridFactory, fileName, verbose, boundarySegments); + + ParentType::gridPtr() = std::shared_ptr<Grid>(gridFactory->createGrid()); + } + } + } +}; + +#endif // HAVE_DUNE_ALUGRID + +#if HAVE_DUNE_FOAMGRID + +/*! + * \brief Provides a grid manager for FoamGrids + * from information in the input file + * + * All keys are expected to be in group GridParameterGroup. + + * The following keys are recognized: + * - File : A DGF or gmsh file to load from, type detection by file extension + * - Verbosity : whether the grid construction should output to standard out + * - LowerLeft : lowerLeft corner of a structured grid + * - UpperRight : upperright corner of a structured grid + * - Cells : number of elements in a structured grid + * + */ +template<int dim, int dimworld> +class GridManager<Dune::FoamGrid<dim, dimworld>> +: public GridManagerBase<Dune::FoamGrid<dim, dimworld>> +{ +public: + using Grid = Dune::FoamGrid<dim, dimworld>; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + // try to create it from file + if (haveParamInGroup(modelParamGroup, "Grid.File")) + { + ParentType::makeGridFromFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"), modelParamGroup); + ParentType::maybeRefineGrid(modelParamGroup); + ParentType::loadBalance(); + return; + } + + // Then look for the necessary keys to construct a structured grid from the input file + else if (haveParamInGroup(modelParamGroup, "Grid.UpperRight")) + { + ParentType::template makeStructuredGrid<dim, dimworld>(ParentType::CellType::Simplex, modelParamGroup); + ParentType::maybeRefineGrid(modelParamGroup); + ParentType::loadBalance(); + } + + // Didn't find a way to construct the grid + else + { + 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"); + + } + } +}; + +/*! + * \brief Provides a grid manager for FoamGrids of dim 1 + * from information in the input file + * + * All keys are expected to be in group GridParameterGroup. + + * The following keys are recognized: + * - File : A DGF or gmsh file to load from, type detection by file extension + * - Verbosity : whether the grid construction should output to standard out + * - LowerLeft : lowerLeft corner of a structured grid + * - UpperRight : upperright corner of a structured grid + * - Cells : number of elements in a structured grid + * + */ +template<int dimworld> +class GridManager<Dune::FoamGrid<1, dimworld>> +: public GridManagerBase<Dune::FoamGrid<1, dimworld>> +{ +public: + using Grid = Dune::FoamGrid<1, dimworld>; + using ParentType = GridManagerBase<Grid>; + + /*! + * \brief Make the grid. This is implemented by specializations of this method. + */ + void init(const std::string& modelParamGroup = "") + { + // try to create it from file + if (haveParamInGroup(modelParamGroup, "Grid.File")) + { + ParentType::makeGridFromFile(getParamFromGroup<std::string>(modelParamGroup, "Grid.File"), modelParamGroup); + ParentType::maybeRefineGrid(modelParamGroup); + ParentType::loadBalance(); + return; + } + + // The required parameters + 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)); + using CellArray = std::array<unsigned int, 1>; + const auto cells = getParamFromGroup<CellArray>(modelParamGroup, "Grid.Cells", std::array<unsigned int, 1>{{1}}); + + // make the grid (structured interval grid in dimworld space) + Dune::GridFactory<Grid> factory; + + constexpr auto geomType = Dune::GeometryTypes::line; + + // create a step vector + GlobalPosition step = upperRight; + step -= lowerLeft, step /= cells[0]; + + // create the vertices + GlobalPosition globalPos = lowerLeft; + for (unsigned int vIdx = 0; vIdx <= cells[0]; vIdx++, globalPos += step) + factory.insertVertex(globalPos); + + // create the cells + for(unsigned int eIdx = 0; eIdx < cells[0]; eIdx++) + factory.insertElement(geomType, {eIdx, eIdx+1}); + + ParentType::gridPtr() = std::shared_ptr<Grid>(factory.createGrid()); + ParentType::maybeRefineGrid(modelParamGroup); + ParentType::loadBalance(); + } +}; + +#endif // HAVE_DUNE_FOAMGRID + +} // end namespace Dumux + +#endif diff --git a/dumux/io/gridcreator.hh b/dumux/io/gridcreator.hh index 36d866273a413cffc753a49a6ab9b19735540c65..daf54134b258833dde2dfccb0d2659be13957201 100644 --- a/dumux/io/gridcreator.hh +++ b/dumux/io/gridcreator.hh @@ -72,68 +72,10 @@ #include <dumux/common/parameters.hh> #include <dumux/discretization/methods.hh> -namespace Dumux -{ +#include "./grid/gmshgriddatahandle.hh" -template<class GridPtr, class GridCreator> -struct GridDataHandle : public Dune::CommDataHandleIF<GridDataHandle<GridPtr, GridCreator>, int> +namespace Dumux { - using GridType = typename GridPtr::element_type; - GridDataHandle( GridPtr& gridPtr) - : gridPtr_(gridPtr), idSet_(gridPtr->localIdSet()) - { - auto gridView = gridPtr_->levelGridView(0); - const auto& indexSet = gridView.indexSet(); - - for( const auto& element : elements(gridView, Dune::Partitions::interior)) - std::swap(GridCreator::elementMarkers_[GridCreator::gridFactory().insertionIndex(element)], data_[idSet_.id(element)]); - - for (const auto& face : entities(gridView, Dune::Codim<1>())) - std::swap(GridCreator::faceMarkers_[indexSet.index(face)], data_[idSet_.id(face)]); - } - - ~GridDataHandle() - { - 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)], data_[idSet_.id(element)]); - - GridCreator::faceMarkers_.resize(indexSet.size(1)); - for (const auto& face : entities(gridView, Dune::Codim<1>())) - std::swap(GridCreator::faceMarkers_[indexSet.index(face)], data_[idSet_.id(face)]); - } - - Dune::CommDataHandleIF<GridDataHandle<GridPtr, GridCreator>, int>& interface() - { return *this; } - - bool contains (int dim, int codim) const - { return codim == 0 || codim == 1; } - - bool fixedSize (int dim, int codim) const - { return true; } - - template<class EntityType> - size_t size (const EntityType& e) const - { return 1; } - - template<class MessageBufferImp, class EntityType> - void gather (MessageBufferImp& buff, const EntityType& e) const - { buff.write(data_[idSet_.id(e)]); } - - template<class MessageBufferImp, class EntityType> - void scatter (MessageBufferImp& buff, const EntityType& e, size_t n) - { buff.read(data_[idSet_.id(e)]); } - -private: - GridPtr &gridPtr_; - using IdSet = typename GridType::LocalIdSet; - const IdSet &idSet_; - mutable std::map< typename IdSet::IdType, int> data_; -}; - /*! * \ingroup InputOutput @@ -145,7 +87,6 @@ class GridCreatorBase { using Intersection = typename Grid::LeafIntersection; using Element = typename Grid::template Codim<0>::Entity; - friend class GridDataHandle<std::shared_ptr<Grid>, GridCreatorBase>; public: /*! * \brief Make the grid. Implement this method in the specialization of this class for a grid type. @@ -315,25 +256,11 @@ public: { { // element and face markers are communicated during load balance - GridDataHandle<std::shared_ptr<Grid>, GridCreatorBase<Grid>> dh(gridPtr()); + using DataHandle = GmshGridDataHandle< Grid, Dune::GridFactory<Grid>, std::vector<int>>; + DataHandle dh(grid(), gridFactory(), elementMarkers_, boundaryMarkers_, faceMarkers_); gridPtr()->loadBalance(dh.interface()); gridPtr()->communicate(dh.interface(), Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication); } - - // transform face markers to boundary markers - boundaryMarkers_.resize(gridPtr()->numBoundarySegments(), 0); - const auto& indexSet = gridPtr()->leafGridView().indexSet(); - for (const auto& element : elements(gridPtr()->leafGridView())) - { - for (const auto& intersection : intersections(gridPtr()->leafGridView(), element)) - { - if (intersection.boundary()) - { - auto marker = faceMarkers_[indexSet.index(element.template subEntity<1>(intersection.indexInInside()))]; - boundaryMarkers_[intersection.boundarySegmentIndex()] = marker; - } - } - } } else gridPtr()->loadBalance(); diff --git a/test/freeflow/navierstokes/test_donea.cc b/test/freeflow/navierstokes/test_donea.cc index 03f60731c4929c6492766f43737d59f6725b6f7f..5e062798d60172861a98d1bb26a149f2f2130fa5 100644 --- a/test/freeflow/navierstokes/test_donea.cc +++ b/test/freeflow/navierstokes/test_donea.cc @@ -52,6 +52,8 @@ #include <dumux/io/staggeredvtkoutputmodule.hh> +#include <dumux/io/grid/gridmanager.hh> + /*! * \brief Provides an interface for customizing error messages associated with * reading in parameters. @@ -102,16 +104,18 @@ int main(int argc, char** argv) try Parameters::init(argc, argv, usage); // try to create a grid (from the given grid file or the input file) - using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); - GridCreator::makeGrid(); - GridCreator::loadBalance(); + using GridManager = Dumux::GridManager<typename GET_PROP_TYPE(TypeTag, Grid), GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>; + GridManager gridMgr; + gridMgr.makeGrid(); + const auto& grid = gridMgr.grid(); + gridMgr.loadBalance(); //////////////////////////////////////////////////////////// // run instationary non-linear problem on this grid //////////////////////////////////////////////////////////// // we compute on the leaf grid view - const auto& leafGridView = GridCreator::grid().leafGridView(); + const auto& leafGridView = grid.leafGridView(); // create the finite volume grid geometry using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); diff --git a/test/io/gridcreator/gridcreatortests.hh b/test/io/gridcreator/gridcreatortests.hh index 8db38bfbc999f15eb892b9b6c67a03c0780c3231..e88caf62a55ffbd8f44bb729aec418eef6c6fa6f 100644 --- a/test/io/gridcreator/gridcreatortests.hh +++ b/test/io/gridcreator/gridcreatortests.hh @@ -30,7 +30,7 @@ #include <dune/grid/common/mcmgmapper.hh> #include <dune/grid/io/file/vtk.hh> -#include <dumux/io/gridcreator.hh> +#include <dumux/io/grid/gridmanager.hh> #include <dumux/discretization/methods.hh> namespace Dumux { @@ -105,7 +105,7 @@ class GridCreatorTests using GridView = typename Grid::LeafGridView; using Scalar = double; static const int dim = Grid::dimension; - using GridCreator = typename Dumux::GridCreatorImpl<Grid, DiscretizationMethod::none>; + using GridManager = typename Dumux::GridManager<Grid>; using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>; public: @@ -114,17 +114,21 @@ public: const std::string& vtkFileName = "test", bool refine = true) { - // initialize the grid - initialize_(); + // make the grid manager and initialize the grid + GridManager gridManager; + gridManager.makeGrid(); + gridManager.loadBalance(); + auto gridData = gridManager.getGridData(); + const auto& leafGridView = gridManager.grid().leafGridView(); // read the boundary and element markers as well as the rank std::vector<int> boundaryMarker, elementMarker, rank; - getBoundaryMarkers_(boundaryMarker); - getElementMarkers_(elementMarker, type); - getRank_(rank); + getBoundaryMarkers_(leafGridView, gridData, boundaryMarker); + getElementMarkers_(leafGridView, gridData, elementMarker, type); + getRank_(leafGridView, rank); // construct a vtk output writer and attach the boundary markers - Dune::VTKSequenceWriter<typename Grid::LeafGridView> vtkWriter(GridCreator::grid().leafGridView(), vtkFileName, ".", ""); + Dune::VTKSequenceWriter<typename Grid::LeafGridView> vtkWriter(leafGridView, vtkFileName, ".", ""); vtkWriter.addVertexData(boundaryMarker, "boundaryMarker"); vtkWriter.addCellData(elementMarker, "elementMarker"); vtkWriter.addCellData(rank, "rank"); @@ -133,10 +137,10 @@ public: if (refine) { // refine grid once and write out the markers again - GridCreator::grid().globalRefine(1); - getBoundaryMarkers_(boundaryMarker); - getElementMarkers_(elementMarker, type); - getRank_(rank); + gridManager.grid().globalRefine(1); + getBoundaryMarkers_(leafGridView, gridData, boundaryMarker); + getElementMarkers_(leafGridView, gridData, elementMarker, type); + getRank_(leafGridView, rank); vtkWriter.write(1); } } @@ -145,16 +149,19 @@ public: const std::string& vtkFileName = "test", bool refine = true) { - // initialize the grid - initialize_(); + // make the grid manager and initialize the grid + GridManager gridManager; + gridManager.makeGrid(); + gridManager.loadBalance(); + auto gridData = gridManager.getGridData(); // read the element markers and the rank std::vector<int> elementMarker, rank; - getElementMarkers_(elementMarker, type); - getRank_(rank); + getElementMarkers_(gridManager.grid().leafGridView(), gridData, elementMarker, type); + getRank_(gridManager.grid().leafGridView(), rank); // construct a vtk output writer and attach the element markers - Dune::VTKSequenceWriter<typename Grid::LeafGridView> vtkWriter(GridCreator::grid().leafGridView(), vtkFileName, ".", ""); + Dune::VTKSequenceWriter<typename Grid::LeafGridView> vtkWriter(gridManager.grid().leafGridView(), vtkFileName, ".", ""); vtkWriter.addCellData(elementMarker, "elementMarker"); vtkWriter.addCellData(rank, "rank"); vtkWriter.write(0); @@ -162,19 +169,17 @@ public: if (refine) { // refine grid once and write out the markers again - GridCreator::grid().globalRefine(1); - getElementMarkers_(elementMarker, type); - getRank_(rank); + gridManager.grid().globalRefine(1); + getElementMarkers_(gridManager.grid().leafGridView(), gridData, elementMarker, type); + getRank_(gridManager.grid().leafGridView(), rank); vtkWriter.write(1); } } private: - static void getRank_(std::vector<int>& rank) + static void getRank_(const GridView& gridView, std::vector<int>& rank) { - const auto& gridView = GridCreator::grid().leafGridView(); - rank.clear(); rank.resize(gridView.size(0)); @@ -185,31 +190,32 @@ private: } } - static void getElementMarkers_(std::vector<int>& elementMarker, + template<class GridData> + static void getElementMarkers_(const GridView& gridView, + const GridData& gridData, + std::vector<int>& elementMarker, const std::string& type) { - const auto& gridView = GridCreator::grid().leafGridView(); - elementMarker.clear(); elementMarker.resize(gridView.size(0)); for(const auto& element : elements(gridView)) { - auto eIdx = gridView.indexSet().index(element); - + const auto eIdx = gridView.indexSet().index(element); if (type == "gmsh") - elementMarker[eIdx] = GridCreator::getElementDomainMarker(element); + elementMarker[eIdx] = gridData->getElementDomainMarker(element); else if (type == "dgf") - elementMarker[eIdx] = GridCreator::parameters(element)[0]; + elementMarker[eIdx] = gridData->parameters(element)[0]; else DUNE_THROW(Dune::InvalidStateException, "No parameters for type " << type); } } - static void getBoundaryMarkers_(std::vector<int>& boundaryMarker) + template<class GridData> + static void getBoundaryMarkers_(const GridView& gridView, + const GridData& gridData, + std::vector<int>& boundaryMarker) { - const auto& gridView = GridCreator::grid().leafGridView(); - boundaryMarker.clear(); boundaryMarker.resize(gridView.size(dim), 0); @@ -217,7 +223,7 @@ private: { for(const auto& intersection : intersections(gridView, element)) { - if(!intersection.boundary()) + if (!intersection.boundary()) continue; const auto refElement = ReferenceElements::general(element.geometry().type()); @@ -231,11 +237,11 @@ private: // make sure we always take the lowest non-zero marker (problem dependent!) if (boundaryMarker[vIdxGlobal] == 0) - boundaryMarker[vIdxGlobal] = GridCreator::getBoundaryDomainMarker(intersection); + boundaryMarker[vIdxGlobal] = gridData->getBoundaryDomainMarker(intersection); else { - if (boundaryMarker[vIdxGlobal] > GridCreator::getBoundaryDomainMarker(intersection)) - boundaryMarker[vIdxGlobal] = GridCreator::getBoundaryDomainMarker(intersection); + if (boundaryMarker[vIdxGlobal] > gridData->getBoundaryDomainMarker(intersection)) + boundaryMarker[vIdxGlobal] = gridData->getBoundaryDomainMarker(intersection); } } } @@ -253,15 +259,6 @@ private: Dune::ForwardCommunication); } } - - static void initialize_() - { - // Make the grid - GridCreator::makeGrid(); - - // Load balancing if parallel - GridCreator::loadBalance(); - } }; } // end namespace Dumux