diff --git a/doc/doxygen/extradoc/featurelist.txt b/doc/doxygen/extradoc/featurelist.txt index 091b69f3641a0879a454807cf87e1685f9644fd3..5030f9041bc9292073da18e36f7c2f3a563b7b32 100644 --- a/doc/doxygen/extradoc/featurelist.txt +++ b/doc/doxygen/extradoc/featurelist.txt @@ -27,7 +27,7 @@ * | :- | :- | * | \c container/test_container_io | Test for writing and reading sequence container to and from file | * | \c gnuplotinterface/test_gnuplotinterface | tests the plotting of data sets, functions, and files | - * | \c gridcreator/test_gridcreator_gmsh | Test for gmsh interface of the grid creator | + * | \c gridmanager/test_gridmanager_gmsh | Test for gmsh interface of the grid manager | * * \section Implicit * | Tests | | Fluidsystem | Gridmanager | Mass/Mole | Non-/Isothermal | Discretization | linearSolver | AdaptiveGrid | BoundaryCondition | Homogeneity | MaterialLaw | Permeability | Gravity | JacobianRecycling | PartialReassemble | Comments | diff --git a/dumux/io/grid/CMakeLists.txt b/dumux/io/grid/CMakeLists.txt index 304ac5b43628f4fa3d4f07965ee915a436f73863..c1118f4cbad1bf7108fccbdd2f7b08ca95151c3a 100644 --- a/dumux/io/grid/CMakeLists.txt +++ b/dumux/io/grid/CMakeLists.txt @@ -1,8 +1,8 @@ install(FILES -cakegridcreator.hh +cakegridmanager.hh cpgridmanager.hh gmshgriddatahandle.hh griddata.hh gridmanager.hh -subgridgridcreator.hh +subgridmanager.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/io/grid) diff --git a/dumux/io/grid/cakegridcreator.hh b/dumux/io/grid/cakegridcreator.hh index 4fe83bc62b7a18b796591f3cd3f314cccedccbbc..b132e9830e48bc421a0b517fa2aa90a1a4b41661 100644 --- a/dumux/io/grid/cakegridcreator.hh +++ b/dumux/io/grid/cakegridcreator.hh @@ -24,15 +24,10 @@ #ifndef DUMUX_CAKE_GRID_CREATOR_HH #define DUMUX_CAKE_GRID_CREATOR_HH -#include <vector> -#include <iostream> -#include <cmath> -#include <algorithm> +#warning "This header is deprecated. Use the new cakegridmanager." -#include <dune/common/dynvector.hh> -#include <dune/common/float_cmp.hh> -#include <dune/grid/common/gridfactory.hh> -#include <dumux/common/parameters.hh> +#include <dune/common/deprecated.hh> +#include "cakegridmanager.hh" namespace Dumux { @@ -42,456 +37,8 @@ namespace Dumux { * with polar Coordinates and one for creating a cartesian grid from * these polar coordinates. */ -template <class Grid> -class CakeGridCreator -{ - using Scalar = typename Grid::ctype; - - using GridFactory = Dune::GridFactory<Grid>; - using GridPointer = std::shared_ptr<Grid>; - - enum { dim = Grid::dimension, - dimWorld = Grid::dimensionworld }; - -public: - /*! - * \brief Make the grid. - */ - void init(const std::string& modelParamGroup = "") - { - static_assert(dim == 2 || dim == 3, "The CakeGridCreator is only implemented for 2D and 3D."); - - const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false); - - std::array<std::vector<Scalar>, dim> polarCoordinates; - // Indices specifing in which direction the piece of cake is oriented - Dune::FieldVector<int, dim> indices(-1); - createVectors(polarCoordinates, indices, modelParamGroup, verbose); - - gridPtr() = createCakeGrid(polarCoordinates, indices, modelParamGroup, verbose); - loadBalance(); - } - - /*! - * \brief Create vectors containing polar coordinates of all points. - * - * All keys are expected to be in group GridParameterGroup. - * The following keys are recognized: - * - Radial : min/max value for radial coordinate - * - Angular : min/max value for angular coordinate - * - Axial : min/max value for axial coordinate - * Adding 0, 1 (or 2 in 3D) specifies in which direction (x, y and z, respectively) - * the radial, angular and axial direction are oriented - * - Cells : number of cells array for x-coordinate (Again, an added 0, 1 or 3 specifies x, y or z - * - Grading : grading factor array for x-coordinate (Same here) - * - Verbosity : whether the grid construction should output to standard out - * - * 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. - */ - static void createVectors(std::array<std::vector<Scalar>, dim> &polarCoordinates, - Dune::FieldVector<int, dim> &indices, - const std::string& modelParamGroup, - bool verbose = false) - { - // The positions - std::array<std::vector<Scalar>, dim> positions; - - for (int i = 0; i < dim; ++i) - { - const bool hasRadial = hasParamInGroup(modelParamGroup, "Grid.Radial" + std::to_string(i)); - const bool hasAngular = hasParamInGroup(modelParamGroup, "Grid.Angular" + std::to_string(i)); - const bool hasAxial = (dim == 3) && hasParamInGroup(modelParamGroup, "Grid.Axial" + std::to_string(i)); - if (static_cast<int>(hasRadial) + static_cast<int>(hasAngular) + static_cast<int>(hasAxial) != 1) - DUNE_THROW(Dune::RangeError, "Multiple or no position vectors (radial, angular, axial) specified in coord direction: " << i << std::endl); - - if (hasRadial) - { - positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Radial" + std::to_string(i)); - indices[0] = i; // Index specifing radial direction - } - - else if (hasAngular) - { - positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Angular" + std::to_string(i)); - indices[1] = i; // Index specifing angular direction - } - - else // hasAxial - { - positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Axial" + std::to_string(i)); - indices[2] = i; // Index specifing axial direction - } - - if (!std::is_sorted(positions[i].begin(), positions[i].end())) - DUNE_THROW(Dune::GridError, "Make sure to specify a monotone increasing \"Positions\" array"); - } - - // check that all indices are specified i.e. > 0 - for (int i = 0; i < dim; ++i) - if (indices[i] < 0) - DUNE_THROW(Dune::RangeError, "Please specify Positions Angular and Radial and Axial correctly and unambiguously!" << std::endl); - - // 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]); - if (cells[i].size() + 1 != positions[i].size()) - DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Cells\" and \"Positions\" arrays"); - } - - // grading factor (has a default) - std::array<std::vector<Scalar>, dim> grading; - for (int i = 0; i < dim; ++i) - { - grading[i].resize(positions[i].size()-1, 1.0); - grading[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]); - if (grading[i].size() + 1 != positions[i].size()) - DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Grading\" and \"Positions\" arrays"); - } - - // make the grid - std::array<std::vector<Scalar>, dim> globalPositions; - using std::pow; - for (int dimIdx = 0; dimIdx < dim; dimIdx++) - { - std::size_t numGlobalPositions = 0; - for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) - numGlobalPositions += cells[dimIdx][zoneIdx] + 1; - - globalPositions[dimIdx].resize(numGlobalPositions); - std::size_t posIdx = 0; - for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) - { - const Scalar lower = positions[dimIdx][zoneIdx]; - const Scalar upper = positions[dimIdx][zoneIdx+1]; - const int numCells = cells[dimIdx][zoneIdx]; - Scalar gradingFactor = grading[dimIdx][zoneIdx]; - const Scalar length = upper - lower; - Scalar 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; - } - } - - // the positions inside a global segment - Dune::DynamicVector<Scalar> localPositions(numCells, 0.0); - for (int i = 1; i < numCells; i++) - { - Scalar hI = height; - if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7)) - { - if (increasingCellSize) - hI *= pow(gradingFactor, i-1); - - else - hI *= pow(gradingFactor, numCells-i); - } - localPositions[i] = localPositions[i-1] + hI; - } - - localPositions *= length; - localPositions += lower; - - for (int i = 0; i < numCells; ++i) - globalPositions[dimIdx][posIdx++] = localPositions[i]; - } - - globalPositions[dimIdx][posIdx] = positions[dimIdx].back(); - } - - polarCoordinates[0] = globalPositions[indices[0]]; - polarCoordinates[1] = globalPositions[indices[1]]; - if (dim == 3) - polarCoordinates[2] = globalPositions[dim - indices[0] - indices[1]]; - - // convert angular coordinates into radians - std::transform(polarCoordinates[1].begin(), polarCoordinates[1].end(), - polarCoordinates[1].begin(), - [](Scalar s){ return s*M_PI/180; }); - } - - /*! - * \brief Creates cartesian grid from polar coordinates. - * - * \param polarCoordinates Vector containing radial, angular and axial coordinates (in this order) - * \param indices Indices specifing the radial, angular and axial direction (in this order) - * \param modelParamGroup name of the model parameter group - * \param verbose if the output should be verbose - */ - std::unique_ptr<Grid> createCakeGrid(std::array<std::vector<Scalar>, dim> &polarCoordinates, - Dune::FieldVector<int, dim> &indices, - const std::string& modelParamGroup, - bool verbose = false) - { - std::vector<Scalar> dR = polarCoordinates[0]; - std::vector<Scalar> dA = polarCoordinates[1]; - - // For the special case of 36o°, the last element is connected with the first one. - // Thus, one needs to distinguish here between all other cases, were the normal procedure - // is applied for all elements until the last one (= dA.size() - 1) and the 360° case, - // where it is only applied until the second to last element (dA.size() - 2). - int maxdA = dA.size() - 1; - if (Dune::FloatCmp::eq(polarCoordinates[1][polarCoordinates[1].size()-1], 360.0)) - { - maxdA = dA.size() - 2; - } - - GridFactory gridFactory; - constexpr auto type = Dune::GeometryTypes::cube(dim); - - // create nodes - if (dim == 3) - { - std::vector<Scalar> dZ = polarCoordinates[2]; - for (int j = 0; j <= dA.size() - 1; ++j) - { - for (int l = 0; l <= dZ.size() - 1; ++l) - { - for (int i = 0; i <= dR.size()- 1; ++i) - { - // Get radius for the well (= a hole) in the center - const auto wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius"); - - // transform into cartesian Coordinates - using std::cos; - using std::sin; - Dune::FieldVector <double, dim> v(0.0); - v[indices[2]] = dZ[l]; - v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i]; - v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i]; - if (verbose) - { - std::cout << "Coordinates of : " - << v[0] << " " << v[1] << " " << v[2] << std::endl; - } - gridFactory.insertVertex(v); - } - } - } - - if (verbose) - std::cout << "Filled node vector" << std::endl; - - // assign nodes - unsigned int z = 0; - unsigned int t = 0; - for (int j = 0; j < dA.size() - 1; ++j) - { - for (int l = 0; l < dZ.size() - 1; ++l) - { - if (j < maxdA) - { - for (int i = 0; i < dR.size() - 1; ++i) - { - unsigned int rSize = dR.size(); - unsigned int zSize = dZ.size(); - std::vector<unsigned int> vid({z, z+1, z+rSize*zSize, - z+rSize*zSize+1, z+rSize, z+rSize+1, - z+rSize*zSize+rSize, z+rSize*zSize+rSize+1}); - - if (verbose) - { - std::cout << "element vertex ids: "; - for (int k = 0; k < vid.size(); ++k) - std::cout << vid[k] << " "; - std::cout << std::endl; - } - - gridFactory.insertElement(type, vid); - - z = z+1; - } - z = z+1; - } - else - { - // assign nodes for 360°-cake - for (int i = 0; i < dR.size() - 1; ++i) - { - // z = z + 1; - unsigned int rSize = dR.size(); - std::vector<unsigned int> vid({z, z+1, t, - t+1, z+rSize, z+rSize+1, - t+rSize, t+rSize+1}); - - if (verbose) - { - std::cout << "element vertex ids: "; - for (int k = 0; k < vid.size(); ++k) - std::cout << vid[k] << " "; - std::cout << std::endl; - } - - gridFactory.insertElement(type, vid); - t = t + 1; - z = z+1; - } - t = t + 1; - z = z+1; - } - if (verbose) - std::cout << "assign nodes 360° ends..." << std::endl; - } - - z = z + dR.size(); - } - } - - // for dim = 2 - else - { - for (int j = 0; j <= dA.size() - 1; ++j) - { - for (int i = 0; i <= dR.size()- 1; ++i) - { - // Get radius for the well (= a hole) in the center - const Scalar wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius"); - - // transform into cartesian Coordinates - Dune::FieldVector <double, dim> v(0.0); - - v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i]; - v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i]; - if(verbose) std::cout << "Coordinates of : " << v[0] << " " << v[1] << std::endl; - gridFactory.insertVertex(v); - } - } - std::cout << "Filled node vector" << std::endl; - - // assign nodes - unsigned int z = 0; - unsigned int t = 0; - for (int j = 0; j < dA.size() - 1; ++j) - { - if (j < maxdA) - { - for (int i = 0; i < dR.size() - 1; ++i) - { - unsigned int rSize = dR.size(); - std::vector<unsigned int> vid({z, z+1, z+rSize, z+rSize+1}); - - if (verbose) - { - std::cout << "element vertex ids: "; - for (int k = 0; k < vid.size(); ++k) - std::cout << vid[k] << " "; - std::cout << std::endl; - } - - gridFactory.insertElement(type, vid); - z = z+1; - } - z = z+1; - } - else - { - // assign nodes for 360°-cake - for (int i = 0; i < dR.size() - 1; ++i) - { - // z = z + 1; - std::vector<unsigned int> vid({z, z+1, t, t+1}); - - if (verbose) - { - std::cout << "element vertex ids: "; - for (int k = 0; k < vid.size(); ++k) - std::cout << vid[k] << " "; - std::cout << std::endl; - } - - gridFactory.insertElement(type, vid); - t = t + 1; - z = z+1; - } - t = t + 1; - z = z+1; - } - std::cout << "assign nodes 360 ends..." << std::endl; - } - } - // return the grid pointer - return std::unique_ptr<Grid>(gridFactory.createGrid()); - } - - /*! - * \brief Returns a reference to the grid. - */ - Grid& grid() - { - return *gridPtr(); - } - - /*! - * \brief Distributes the grid on all processes of a parallel - * computation. - */ - void loadBalance() - { - gridPtr()->loadBalance(); - } - -protected: - - /*! - * \brief Returns a reference to the shared pointer to the grid. - */ - GridPointer& gridPtr() - { - return cakeGrid_; - } - -private: - GridPointer cakeGrid_; -}; - +template<class Grid> +using CakeGridCreator DUNE_DEPRECATED_MSG("Use CakeGridManager instead!") = CakeGridManager<Grid>; } // end namespace Dumux #endif diff --git a/dumux/io/grid/cakegridmanager.hh b/dumux/io/grid/cakegridmanager.hh new file mode 100644 index 0000000000000000000000000000000000000000..0e64861954c7b5988a54dd66f365f1ae9c6878e8 --- /dev/null +++ b/dumux/io/grid/cakegridmanager.hh @@ -0,0 +1,497 @@ +// -*- 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 3 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 a piece of cake grid + */ +#ifndef DUMUX_CAKE_GRID_MANAGER_HH +#define DUMUX_CAKE_GRID_MANAGER_HH + +#include <vector> +#include <iostream> +#include <cmath> +#include <algorithm> + +#include <dune/common/dynvector.hh> +#include <dune/common/float_cmp.hh> +#include <dune/grid/common/gridfactory.hh> +#include <dumux/common/parameters.hh> + +namespace Dumux { + +/*! + * \ingroup InputOutput + * \brief Provides a grid manager with a method for creating creating vectors + * with polar Coordinates and one for creating a cartesian grid from + * these polar coordinates. + */ +template <class Grid> +class CakeGridManager +{ + using Scalar = typename Grid::ctype; + + using GridFactory = Dune::GridFactory<Grid>; + using GridPointer = std::shared_ptr<Grid>; + + enum { dim = Grid::dimension, + dimWorld = Grid::dimensionworld }; + +public: + /*! + * \brief Make the grid. + */ + void init(const std::string& modelParamGroup = "") + { + static_assert(dim == 2 || dim == 3, "The CakeGridManager is only implemented for 2D and 3D."); + + const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false); + + std::array<std::vector<Scalar>, dim> polarCoordinates; + // Indices specifing in which direction the piece of cake is oriented + Dune::FieldVector<int, dim> indices(-1); + createVectors(polarCoordinates, indices, modelParamGroup, verbose); + + gridPtr() = createCakeGrid(polarCoordinates, indices, modelParamGroup, verbose); + loadBalance(); + } + + /*! + * \brief Create vectors containing polar coordinates of all points. + * + * All keys are expected to be in group GridParameterGroup. + * The following keys are recognized: + * - Radial : min/max value for radial coordinate + * - Angular : min/max value for angular coordinate + * - Axial : min/max value for axial coordinate + * Adding 0, 1 (or 2 in 3D) specifies in which direction (x, y and z, respectively) + * the radial, angular and axial direction are oriented + * - Cells : number of cells array for x-coordinate (Again, an added 0, 1 or 3 specifies x, y or z + * - Grading : grading factor array for x-coordinate (Same here) + * - Verbosity : whether the grid construction should output to standard out + * + * 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. + */ + static void createVectors(std::array<std::vector<Scalar>, dim> &polarCoordinates, + Dune::FieldVector<int, dim> &indices, + const std::string& modelParamGroup, + bool verbose = false) + { + // The positions + std::array<std::vector<Scalar>, dim> positions; + + for (int i = 0; i < dim; ++i) + { + const bool hasRadial = hasParamInGroup(modelParamGroup, "Grid.Radial" + std::to_string(i)); + const bool hasAngular = hasParamInGroup(modelParamGroup, "Grid.Angular" + std::to_string(i)); + const bool hasAxial = (dim == 3) && hasParamInGroup(modelParamGroup, "Grid.Axial" + std::to_string(i)); + if (static_cast<int>(hasRadial) + static_cast<int>(hasAngular) + static_cast<int>(hasAxial) != 1) + DUNE_THROW(Dune::RangeError, "Multiple or no position vectors (radial, angular, axial) specified in coord direction: " << i << std::endl); + + if (hasRadial) + { + positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Radial" + std::to_string(i)); + indices[0] = i; // Index specifing radial direction + } + + else if (hasAngular) + { + positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Angular" + std::to_string(i)); + indices[1] = i; // Index specifing angular direction + } + + else // hasAxial + { + positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Axial" + std::to_string(i)); + indices[2] = i; // Index specifing axial direction + } + + if (!std::is_sorted(positions[i].begin(), positions[i].end())) + DUNE_THROW(Dune::GridError, "Make sure to specify a monotone increasing \"Positions\" array"); + } + + // check that all indices are specified i.e. > 0 + for (int i = 0; i < dim; ++i) + if (indices[i] < 0) + DUNE_THROW(Dune::RangeError, "Please specify Positions Angular and Radial and Axial correctly and unambiguously!" << std::endl); + + // 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]); + if (cells[i].size() + 1 != positions[i].size()) + DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Cells\" and \"Positions\" arrays"); + } + + // grading factor (has a default) + std::array<std::vector<Scalar>, dim> grading; + for (int i = 0; i < dim; ++i) + { + grading[i].resize(positions[i].size()-1, 1.0); + grading[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]); + if (grading[i].size() + 1 != positions[i].size()) + DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Grading\" and \"Positions\" arrays"); + } + + // make the grid + std::array<std::vector<Scalar>, dim> globalPositions; + using std::pow; + for (int dimIdx = 0; dimIdx < dim; dimIdx++) + { + std::size_t numGlobalPositions = 0; + for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) + numGlobalPositions += cells[dimIdx][zoneIdx] + 1; + + globalPositions[dimIdx].resize(numGlobalPositions); + std::size_t posIdx = 0; + for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) + { + const Scalar lower = positions[dimIdx][zoneIdx]; + const Scalar upper = positions[dimIdx][zoneIdx+1]; + const int numCells = cells[dimIdx][zoneIdx]; + Scalar gradingFactor = grading[dimIdx][zoneIdx]; + const Scalar length = upper - lower; + Scalar 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; + } + } + + // the positions inside a global segment + Dune::DynamicVector<Scalar> localPositions(numCells, 0.0); + for (int i = 1; i < numCells; i++) + { + Scalar hI = height; + if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7)) + { + if (increasingCellSize) + hI *= pow(gradingFactor, i-1); + + else + hI *= pow(gradingFactor, numCells-i); + } + localPositions[i] = localPositions[i-1] + hI; + } + + localPositions *= length; + localPositions += lower; + + for (int i = 0; i < numCells; ++i) + globalPositions[dimIdx][posIdx++] = localPositions[i]; + } + + globalPositions[dimIdx][posIdx] = positions[dimIdx].back(); + } + + polarCoordinates[0] = globalPositions[indices[0]]; + polarCoordinates[1] = globalPositions[indices[1]]; + if (dim == 3) + polarCoordinates[2] = globalPositions[dim - indices[0] - indices[1]]; + + // convert angular coordinates into radians + std::transform(polarCoordinates[1].begin(), polarCoordinates[1].end(), + polarCoordinates[1].begin(), + [](Scalar s){ return s*M_PI/180; }); + } + + /*! + * \brief Creates cartesian grid from polar coordinates. + * + * \param polarCoordinates Vector containing radial, angular and axial coordinates (in this order) + * \param indices Indices specifing the radial, angular and axial direction (in this order) + * \param modelParamGroup name of the model parameter group + * \param verbose if the output should be verbose + */ + std::unique_ptr<Grid> createCakeGrid(std::array<std::vector<Scalar>, dim> &polarCoordinates, + Dune::FieldVector<int, dim> &indices, + const std::string& modelParamGroup, + bool verbose = false) + { + std::vector<Scalar> dR = polarCoordinates[0]; + std::vector<Scalar> dA = polarCoordinates[1]; + + // For the special case of 36o°, the last element is connected with the first one. + // Thus, one needs to distinguish here between all other cases, were the normal procedure + // is applied for all elements until the last one (= dA.size() - 1) and the 360° case, + // where it is only applied until the second to last element (dA.size() - 2). + int maxdA = dA.size() - 1; + if (Dune::FloatCmp::eq(polarCoordinates[1][polarCoordinates[1].size()-1], 360.0)) + { + maxdA = dA.size() - 2; + } + + GridFactory gridFactory; + constexpr auto type = Dune::GeometryTypes::cube(dim); + + // create nodes + if (dim == 3) + { + std::vector<Scalar> dZ = polarCoordinates[2]; + for (int j = 0; j <= dA.size() - 1; ++j) + { + for (int l = 0; l <= dZ.size() - 1; ++l) + { + for (int i = 0; i <= dR.size()- 1; ++i) + { + // Get radius for the well (= a hole) in the center + const auto wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius"); + + // transform into cartesian Coordinates + using std::cos; + using std::sin; + Dune::FieldVector <double, dim> v(0.0); + v[indices[2]] = dZ[l]; + v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i]; + v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i]; + if (verbose) + { + std::cout << "Coordinates of : " + << v[0] << " " << v[1] << " " << v[2] << std::endl; + } + gridFactory.insertVertex(v); + } + } + } + + if (verbose) + std::cout << "Filled node vector" << std::endl; + + // assign nodes + unsigned int z = 0; + unsigned int t = 0; + for (int j = 0; j < dA.size() - 1; ++j) + { + for (int l = 0; l < dZ.size() - 1; ++l) + { + if (j < maxdA) + { + for (int i = 0; i < dR.size() - 1; ++i) + { + unsigned int rSize = dR.size(); + unsigned int zSize = dZ.size(); + std::vector<unsigned int> vid({z, z+1, z+rSize*zSize, + z+rSize*zSize+1, z+rSize, z+rSize+1, + z+rSize*zSize+rSize, z+rSize*zSize+rSize+1}); + + if (verbose) + { + std::cout << "element vertex ids: "; + for (int k = 0; k < vid.size(); ++k) + std::cout << vid[k] << " "; + std::cout << std::endl; + } + + gridFactory.insertElement(type, vid); + + z = z+1; + } + z = z+1; + } + else + { + // assign nodes for 360°-cake + for (int i = 0; i < dR.size() - 1; ++i) + { + // z = z + 1; + unsigned int rSize = dR.size(); + std::vector<unsigned int> vid({z, z+1, t, + t+1, z+rSize, z+rSize+1, + t+rSize, t+rSize+1}); + + if (verbose) + { + std::cout << "element vertex ids: "; + for (int k = 0; k < vid.size(); ++k) + std::cout << vid[k] << " "; + std::cout << std::endl; + } + + gridFactory.insertElement(type, vid); + t = t + 1; + z = z+1; + } + t = t + 1; + z = z+1; + } + if (verbose) + std::cout << "assign nodes 360° ends..." << std::endl; + } + + z = z + dR.size(); + } + } + + // for dim = 2 + else + { + for (int j = 0; j <= dA.size() - 1; ++j) + { + for (int i = 0; i <= dR.size()- 1; ++i) + { + // Get radius for the well (= a hole) in the center + const Scalar wellRadius = getParamFromGroup<Scalar>(modelParamGroup, "Grid.WellRadius"); + + // transform into cartesian Coordinates + Dune::FieldVector <double, dim> v(0.0); + + v[indices[0]] = cos(dA[j])*wellRadius + cos(dA[j])*dR[i]; + v[indices[1]] = sin(dA[j])*wellRadius + sin(dA[j])*dR[i]; + if(verbose) std::cout << "Coordinates of : " << v[0] << " " << v[1] << std::endl; + gridFactory.insertVertex(v); + } + } + std::cout << "Filled node vector" << std::endl; + + // assign nodes + unsigned int z = 0; + unsigned int t = 0; + for (int j = 0; j < dA.size() - 1; ++j) + { + if (j < maxdA) + { + for (int i = 0; i < dR.size() - 1; ++i) + { + unsigned int rSize = dR.size(); + std::vector<unsigned int> vid({z, z+1, z+rSize, z+rSize+1}); + + if (verbose) + { + std::cout << "element vertex ids: "; + for (int k = 0; k < vid.size(); ++k) + std::cout << vid[k] << " "; + std::cout << std::endl; + } + + gridFactory.insertElement(type, vid); + z = z+1; + } + z = z+1; + } + else + { + // assign nodes for 360°-cake + for (int i = 0; i < dR.size() - 1; ++i) + { + // z = z + 1; + std::vector<unsigned int> vid({z, z+1, t, t+1}); + + if (verbose) + { + std::cout << "element vertex ids: "; + for (int k = 0; k < vid.size(); ++k) + std::cout << vid[k] << " "; + std::cout << std::endl; + } + + gridFactory.insertElement(type, vid); + t = t + 1; + z = z+1; + } + t = t + 1; + z = z+1; + } + std::cout << "assign nodes 360 ends..." << std::endl; + } + } + // return the grid pointer + return std::unique_ptr<Grid>(gridFactory.createGrid()); + } + + /*! + * \brief Returns a reference to the grid. + */ + Grid& grid() + { + return *gridPtr(); + } + + /*! + * \brief Distributes the grid on all processes of a parallel + * computation. + */ + void loadBalance() + { + gridPtr()->loadBalance(); + } + +protected: + + /*! + * \brief Returns a reference to the shared pointer to the grid. + */ + GridPointer& gridPtr() + { + return cakeGrid_; + } + +private: + GridPointer cakeGrid_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/io/grid/subgridgridcreator.hh b/dumux/io/grid/subgridgridcreator.hh index 139f2006224dbcae3a41a67447b910ab5a582670..fbe26e1bb30fe98d48b78018dd8d00bb96cccc5a 100644 --- a/dumux/io/grid/subgridgridcreator.hh +++ b/dumux/io/grid/subgridgridcreator.hh @@ -26,13 +26,10 @@ #if HAVE_DUNE_SUBGRID -#include <memory> +#warning "This header is deprecated. Use the new subgridmanager." -#include <dune/subgrid/subgrid.hh> -#include <dune/grid/io/file/vtk.hh> -#include <dune/grid/io/file/dgfparser/dgfwriter.hh> - -#include <dumux/common/parameters.hh> +#include <dune/common/deprecated.hh> +#include "subgridgridmanager.hh" namespace Dumux { @@ -41,63 +38,7 @@ namespace Dumux { * \brief A grid creator for dune-subgrid. */ template <class HostGrid> -class SubgridGridCreator -{ - static constexpr auto dim = HostGrid::dimension; - -public: - using Grid = Dune::SubGrid<dim, HostGrid>; - - /*! - * \brief Make the subgrid. - */ - template<class ElementSelector> - static std::unique_ptr<Grid> makeGrid(HostGrid& hostgrid, - const ElementSelector& selector, - const std::string& modelParamGroup = "") - { - // A unique pointer to the subgrid. - auto subgridPtr = std::make_unique<Grid>(hostgrid); - - // A container to store the host grid elements' ids. - std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid; - const auto& globalIDset = subgridPtr->getHostGrid().globalIdSet(); - - // Construct the subgrid. - subgridPtr->createBegin(); - - // Loop over all elements of the host grid and use the selector to - // choose which elements to add to the subgrid. - auto hostGridView = subgridPtr->getHostGrid().leafGridView(); - for (const auto& e : elements(hostGridView)) - if(selector(e)) - elementsForSubgrid.insert(globalIDset.template id<0>(e)); - - subgridPtr->insertSetPartial(elementsForSubgrid); - subgridPtr->createEnd(); - - // If desired, write out the final subgrid as a dgf file. - if(getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToDGF", false)) - { - const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", ""); - const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix; - Dune::DGFWriter<typename Grid::LeafGridView> writer(subgridPtr->leafGridView()); - writer.write(name + ".dgf"); - } - - // If desired, write out the hostgrid as vtk file. - if(getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToVtk", false)) - { - const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", ""); - const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix; - Dune::VTKWriter<typename Grid::LeafGridView> vtkWriter(subgridPtr->leafGridView()); - vtkWriter.write(name); - } - - // Return a unique pointer to the subgrid. - return subgridPtr; - } -}; +using SubgridGridCreator DUNE_DEPRECATED_MSG("Use SubgridManager instead!") = SubgridManager<HostGrid>; } // end namespace Dumux diff --git a/dumux/io/grid/subgridmanager.hh b/dumux/io/grid/subgridmanager.hh new file mode 100644 index 0000000000000000000000000000000000000000..d87d4b1506ee68159dea481c72da5ac3ad940a8f --- /dev/null +++ b/dumux/io/grid/subgridmanager.hh @@ -0,0 +1,105 @@ +// -*- 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 3 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 grid manager for dune-subgrid. + */ +#ifndef DUMUX_SUBGRID_MANAGER_HH +#define DUMUX_SUBGRID_MANAGER_HH + +#if HAVE_DUNE_SUBGRID + +#include <memory> + +#include <dune/subgrid/subgrid.hh> +#include <dune/grid/io/file/vtk.hh> +#include <dune/grid/io/file/dgfparser/dgfwriter.hh> + +#include <dumux/common/parameters.hh> + +namespace Dumux { + +/*! + * \ingroup InputOutput + * \brief A grid manager for dune-subgrid. + */ +template <class HostGrid> +class SubgridManager +{ + static constexpr auto dim = HostGrid::dimension; + +public: + using Grid = Dune::SubGrid<dim, HostGrid>; + + /*! + * \brief Make the subgrid. + */ + template<class ElementSelector> + static std::unique_ptr<Grid> makeGrid(HostGrid& hostgrid, + const ElementSelector& selector, + const std::string& modelParamGroup = "") + { + // A unique pointer to the subgrid. + auto subgridPtr = std::make_unique<Grid>(hostgrid); + + // A container to store the host grid elements' ids. + std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid; + const auto& globalIDset = subgridPtr->getHostGrid().globalIdSet(); + + // Construct the subgrid. + subgridPtr->createBegin(); + + // Loop over all elements of the host grid and use the selector to + // choose which elements to add to the subgrid. + auto hostGridView = subgridPtr->getHostGrid().leafGridView(); + for (const auto& e : elements(hostGridView)) + if(selector(e)) + elementsForSubgrid.insert(globalIDset.template id<0>(e)); + + subgridPtr->insertSetPartial(elementsForSubgrid); + subgridPtr->createEnd(); + + // If desired, write out the final subgrid as a dgf file. + if(getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToDGF", false)) + { + const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", ""); + const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix; + Dune::DGFWriter<typename Grid::LeafGridView> writer(subgridPtr->leafGridView()); + writer.write(name + ".dgf"); + } + + // If desired, write out the hostgrid as vtk file. + if(getParamFromGroup<bool>(modelParamGroup, "Grid.WriteSubGridToVtk", false)) + { + const auto postfix = getParamFromGroup<std::string>(modelParamGroup, "Problem.Name", ""); + const std::string name = postfix == "" ? "subgrid" : "subgrid_" + postfix; + Dune::VTKWriter<typename Grid::LeafGridView> vtkWriter(subgridPtr->leafGridView()); + vtkWriter.write(name); + } + + // Return a unique pointer to the subgrid. + return subgridPtr; + } +}; + +} // end namespace Dumux + +#endif // HAVE_DUNE_SUBGRID +#endif // DUMUX_SUBGRID_MANAGER_HH diff --git a/test/io/gridmanager/test_gridmanager_cake.cc b/test/io/gridmanager/test_gridmanager_cake.cc index da1c37d3f73a5afb0961ce38513be05d44e462c9..0e202f0ac10dd95349eb558e0788263a23ecc2ca 100644 --- a/test/io/gridmanager/test_gridmanager_cake.cc +++ b/test/io/gridmanager/test_gridmanager_cake.cc @@ -17,7 +17,7 @@ /*! * \file * - * \brief Test for the cake grid creator + * \brief Test for the cake grid manager */ #include<string> @@ -30,7 +30,7 @@ #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> #include <dumux/io/grid/gridmanager.hh> -#include <dumux/io/grid/cakegridcreator.hh> +#include <dumux/io/grid/cakegridmanager.hh> #if HAVE_UG #include <dune/grid/uggrid.hh> @@ -52,10 +52,10 @@ using Grid = Dune::UGGrid<dim>; #endif template<int dim> -void testCakeGridCreator(const std::string& name) +void testCakeGridManager(const std::string& name) { // using declarations - using GridManager = typename Dumux::CakeGridCreator<Grid<dim>>; + using GridManager = typename Dumux::CakeGridManager<Grid<dim>>; GridManager gridManager; // make the grid @@ -78,10 +78,10 @@ int main(int argc, char** argv) try const auto name = Dumux::getParam<std::string>("Grid.Name"); // test 3-D - testCakeGridCreator<3>("cake-3d-" + name); + testCakeGridManager<3>("cake-3d-" + name); // test 2-D - testCakeGridCreator<2>("cake-2d-" + name); + testCakeGridManager<2>("cake-2d-" + name); return 0; } diff --git a/test/io/gridmanager/test_gridmanager_subgrid.cc b/test/io/gridmanager/test_gridmanager_subgrid.cc index 9db016d7019eb9b4132d6d71748e19ce3fdfee55..0884ee3da7d6b9469abbde60b638586498a04044 100644 --- a/test/io/gridmanager/test_gridmanager_subgrid.cc +++ b/test/io/gridmanager/test_gridmanager_subgrid.cc @@ -16,7 +16,7 @@ *****************************************************************************/ /*! * \file - * \brief Test for the cake grid creator + * \brief Test for the cake grid manager */ #include <config.h> #include <iostream> @@ -29,7 +29,7 @@ #include <dumux/common/parameters.hh> #include <dumux/io/grid/gridmanager.hh> -#include <dumux/io/grid/subgridgridcreator.hh> +#include <dumux/io/grid/subgridmanager.hh> #include <dumux/discretization/method.hh> /*! @@ -110,9 +110,9 @@ int main(int argc, char** argv) try CircleSelector<GlobalPosition> elementSelectorThree(center); // Create three different subgrids from the same hostgrid. - auto subgridPtrOne = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorOne, "SubGridOne"); - auto subgridPtrTwo = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorTwo, "SubGridTwo"); - auto subgridPtrThree = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorThree, "SubGridThree"); + auto subgridPtrOne = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorOne, "SubGridOne"); + auto subgridPtrTwo = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorTwo, "SubGridTwo"); + auto subgridPtrThree = SubgridManager<HostGrid>::makeGrid(hostGrid, elementSelectorThree, "SubGridThree"); std::cout << "Constructing a host grid and three subgrids took " << timer.elapsed() << " seconds.\n"; diff --git a/test/multidomain/boundary/darcydarcy/1p_1p/main.cc b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc index 7efd65b12e6bbfb14b905a77ba38e159b2ee151d..790e178a6e578b525fe454449dfa6ff161d81ce1 100644 --- a/test/multidomain/boundary/darcydarcy/1p_1p/main.cc +++ b/test/multidomain/boundary/darcydarcy/1p_1p/main.cc @@ -39,7 +39,7 @@ #include <dumux/discretization/method.hh> #include <dumux/io/vtkoutputmodule.hh> #include <dumux/io/grid/gridmanager.hh> -#include <dumux/io/grid/subgridgridcreator.hh> +#include <dumux/io/grid/subgridmanager.hh> #include <dumux/porousmediumflow/1p/model.hh> #include <dumux/material/components/simpleh2o.hh> @@ -156,8 +156,8 @@ int main(int argc, char** argv) try auto elementSelector1 = [&lensLowerLeft, &lensUpperRight](const auto& element) { return !LensSpatialParams::pointInLens(element.geometry().center(), lensLowerLeft, lensUpperRight); }; - auto subGrid0 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector0); - auto subGrid1 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector1); + auto subGrid0 = SubgridManager<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector0); + auto subGrid1 = SubgridManager<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector1); //////////////////////////////////////////////////////////// // run instationary non-linear problem on this grid diff --git a/test/multidomain/boundary/darcydarcy/1p_2p/main.cc b/test/multidomain/boundary/darcydarcy/1p_2p/main.cc index 623d02f1fa26114cc3530198d628c71c9e71cc85..af9c56bff21044dcff62c172bf7a1d4a217db3ab 100644 --- a/test/multidomain/boundary/darcydarcy/1p_2p/main.cc +++ b/test/multidomain/boundary/darcydarcy/1p_2p/main.cc @@ -39,7 +39,7 @@ #include <dumux/discretization/method.hh> #include <dumux/io/vtkoutputmodule.hh> #include <dumux/io/grid/gridmanager.hh> -#include <dumux/io/grid/subgridgridcreator.hh> +#include <dumux/io/grid/subgridmanager.hh> #include <dumux/porousmediumflow/1p/model.hh> #include <dumux/porousmediumflow/2p/model.hh> @@ -158,8 +158,8 @@ int main(int argc, char** argv) try auto elementSelector1 = [radius](const auto& element) { return element.geometry().center().two_norm() < radius; }; - auto subGrid0 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector0); - auto subGrid1 = SubgridGridCreator<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector1); + auto subGrid0 = SubgridManager<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector0); + auto subGrid1 = SubgridManager<FullDomainGrid>::makeGrid(gridManager.grid(), elementSelector1); //////////////////////////////////////////////////////////// // run instationary non-linear problem on this grid