GraphOfGrid.hpp

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright 2024 Equinor ASA.

  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 3 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 OPM_GRAPH_OF_GRID_HEADER
#define OPM_GRAPH_OF_GRID_HEADER

#include <opm/grid/CpGrid.hpp>

namespace Opm {

template<typename Grid>
class GraphOfGrid{
    using WeightType = float;
    using EdgeList = std::map<int,WeightType>;

    struct VertexProperties
    {
        int nproc = 0; // number of processor
        WeightType weight = 1; // vertex weight
        EdgeList edges;
    };

public:
    explicit GraphOfGrid (const Grid& grid_,
                          const double* transmissibilities=nullptr,
                          const Dune::EdgeWeightMethod edgeWeightMethod=Dune::EdgeWeightMethod::defaultTransEdgeWgt,
                          int level = -1)
        : grid(grid_)
    {
        createGraph(transmissibilities,edgeWeightMethod, level);
    }

    const Grid& getGrid() const
    {
        return grid;
    }

    int size () const
    {
        return graph.size();
    }

    auto begin() const
    {
        return graph.begin();
    }
    auto end() const
    {
        return graph.end();
    }

    auto find(int gID) const
    {
        // search the graph first, and then wells
        auto pgID = graph.find(gID);
        if (pgID == graph.end())
        {
            gID = wellID(gID);
            if (gID == -1)
            {
                return graph.end();
            }
            pgID = graph.find(gID);
        }
        return pgID;
    }

    const VertexProperties& getVertex (int gID) const
    {
        auto pgID = find(gID);
        if (pgID == graph.end())
        {
            OPM_THROW(std::logic_error, "GraphOfGrid::getVertex: gID is not in the graph!");
        }
        return pgID->second;
    }

    // returns -1 if vertex with such global ID is not in the graph (or wells)
    int numEdges (int gID) const
    {
        auto pgID = find(gID);
        if (pgID == graph.end())
        {
            return -1;
        }
        else
        {
            return pgID->second.edges.size();
        }
    }

    const EdgeList& edgeList(int gID) const
    {
        // get iterator to the vertex or the well containing it
        auto pgID = find(gID);
        if (pgID==graph.end())
        {
            OPM_THROW(std::logic_error, "GraphOfGrid::edgeList: gID is not in the graph!");
        }
        return pgID->second.edges;
    }

    int contractVertices (int gID1, int gID2);

    void addWell (const std::set<int>& well, bool checkIntersection=true);

    const auto& getWells () const
    {
        return wells;
    }

    void addNeighboringCellsToWells ();
    void addNeighboringCellsToWells (int layers)
    {
        for (int i=0; i<layers; ++i)
        {
            addNeighboringCellsToWells();
        }
    }

private:
    void createGraph (const double* transmissibilities=nullptr,
                      const Dune::EdgeWeightMethod edgeWeightMethod=Dune::EdgeWeightMethod::defaultTransEdgeWgt,
                      int level = -1);

    int wellID (int gID) const;

    void mergeWellIndices(const std::set<int>& well);

    void contractWellAndAdd(const std::set<int>& well);

    const Grid& grid;
    std::unordered_map<int, VertexProperties> graph; // <gID, VertexProperties>
    std::list<std::set<int>> wells;
};

} // namespace Opm

#endif // OPM_GRAPH_OF_GRID_HEADER