dgfvanguard.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).
 
  OPM 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.
 
  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_DGF_GRID_VANGUARD_HH
#define EWOMS_DGF_GRID_VANGUARD_HH
 
#include <dune/grid/io/file/dgfparser/dgfparser.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include <opm/models/discretefracture/fracturemapper.hh>
#include <opm/models/io/basevanguard.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/parametersystem.hpp>
 
#include <memory>
#include <string>
#include <vector>
 
namespace Opm {
 
template <class TypeTag>
class DgfVanguard : public BaseVanguard<TypeTag>
{
    using ParentType = BaseVanguard<TypeTag>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using FractureMapper = Opm::FractureMapper<TypeTag>;
 
    using GridPointer = std::unique_ptr<Grid>;
 
public:
    static void registerParameters()
    {
        Parameters::Register<Parameters::GridFile>
            ("The file name of the DGF file to load");
        Parameters::Register<Parameters::GridGlobalRefinements>
            ("The number of global refinements of the grid "
             "executed after it was loaded");
    }
 
    explicit DgfVanguard(Simulator& simulator)
        : ParentType(simulator)
    {
        const std::string dgfFileName = Parameters::Get<Parameters::GridFile>();
        const unsigned numRefinments = Parameters::Get<Parameters::GridGlobalRefinements>();
 
        {
            // create DGF GridPtr from a dgf file
            Dune::GridPtr<Grid> dgfPointer(dgfFileName);
 
            // this is only implemented for 2d currently
            addFractures_(dgfPointer);
 
            // store pointer to dune grid
            gridPtr_.reset(dgfPointer.release());
        }
 
        if (numRefinments > 0) {
            gridPtr_->globalRefine(static_cast<int>(numRefinments));
        }
 
        this->finalizeInit_();
    }
 
    Grid& grid()
    { return *gridPtr_; }
 
    const Grid& grid() const
    { return *gridPtr_; }
 
    void loadBalance()
    { gridPtr_->loadBalance(); }
 
    FractureMapper& fractureMapper()
    { return fractureMapper_; }
 
    const FractureMapper& fractureMapper() const
    { return fractureMapper_; }
 
protected:
    void addFractures_(Dune::GridPtr<Grid>& dgfPointer)
    {
        using LevelGridView = typename Grid::LevelGridView;
 
        // check if fractures are available (only 2d currently)
        if (dgfPointer.nofParameters(static_cast<int>(Grid::dimension)) == 0) {
            return;
        }
 
        LevelGridView gridView = dgfPointer->levelGridView(/*level=*/0);
        const unsigned edgeCodim = Grid::dimension - 1;
 
        using VertexMapper = Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView>;
        VertexMapper vertexMapper(gridView, Dune::mcmgVertexLayout());
 
        // first create a map of the dune to ART vertex indices
        for (const auto& element: elements(gridView)) {
            const auto& refElem =
                Dune::ReferenceElements<Scalar, Grid::dimension>::general(element.type());
 
            const int edges = refElem.size(edgeCodim);
            for (int edge = 0; edge < edges; ++edge) {
                const int vertices = refElem.size(edge, edgeCodim, Grid::dimension);
                std::vector<unsigned> vertexIndices;
                vertexIndices.reserve(Grid::dimension);
                for (int vx = 0; vx < vertices; ++vx) {
                    // get local vertex number from edge
                    const int localVx = refElem.subEntity(edge, edgeCodim, vx, Grid::dimension);
 
                    // get vertex
                    const auto vertex = element.template subEntity<Grid::dimension>(localVx);
 
                    // if vertex has parameter 1 insert as a fracture vertex
                    if (dgfPointer.parameters(vertex)[0] > 0) {
                        vertexIndices.push_back(
                            static_cast<unsigned>(vertexMapper.subIndex(element,
                                                                        static_cast<int>(localVx),
                                                                        Grid::dimension)));
                    }
                }
                // if 2 vertices have been found with flag 1 insert a fracture edge
                if (static_cast<int>(vertexIndices.size()) == Grid::dimension) {
                    fractureMapper_.addFractureEdge(vertexIndices[0], vertexIndices[1]);
                }
            }
        }
    }
 
private:
    GridPointer    gridPtr_;
    FractureMapper fractureMapper_;
};
 
} // namespace Opm
 
#endif