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