OrientedEntityTable.hpp

//===========================================================================
//
// File: OrientedEntityTable.hpp
//
// Created: Wed Aug 26 11:13:20 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Bård Skaflestad     <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil 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/>.
*/

#ifndef OPM_ORIENTEDENTITYTABLE_HEADER
#define OPM_ORIENTEDENTITYTABLE_HEADER

#include "EntityRep.hpp"
#include <opm/grid/utility/SparseTable.hpp>
#include <map>
#include <climits>

namespace Dune
{
    namespace cpgrid
    {


        template <int codim_to>
        class OrientedEntityRange : private Opm::SparseTable< EntityRep<codim_to> >::row_type
        {
        public:
            typedef EntityRep<codim_to> ToType;
            typedef ToType* ToTypePtr;
            typedef typename Opm::SparseTable<ToType>::row_type R;

            OrientedEntityRange()
                : orientation_(true)
            {
            }
            OrientedEntityRange(const R& r, bool orientation)
                : R(r), orientation_(orientation)
            {
            }
            int size () const { return R::size(); }
            using R::empty;
            using R::begin;
            using R::end;
            ToType operator[](int subindex) const
            {
                ToType erep = *(this->begin() + subindex);
                return orientation_ ? erep : erep.opposite();
            }
        private:
            bool orientation_;
        };


        template <int codim_to>
        class MutableOrientedEntityRange : private Opm::SparseTable< EntityRep<codim_to> >::mutable_row_type
        {
        public:
            typedef EntityRep<codim_to> ToType;
            typedef ToType* ToTypePtr;
            typedef typename Opm::SparseTable<ToType>::mutable_row_type R;

            MutableOrientedEntityRange()
                : R(ToTypePtr(0), ToTypePtr(0)), orientation_(true)
            {
            }
            MutableOrientedEntityRange(const R& r, bool orientation)
                : R(r), orientation_(orientation)
            {
            }
            int size () const { return R::size(); }
            using R::empty;
            using R::begin;
            using R::end;
            ToType operator[](int subindex) const
            {
                ToType erep = R::operator[](subindex);
                return orientation_ ? erep : erep.opposite();
            }
        private:
            bool orientation_;
        };


        template <int codim_from, int codim_to>
        class OrientedEntityTable : private Opm::SparseTable< EntityRep<codim_to> >
        {
            friend class CpGridData;
        public:
            typedef EntityRep<codim_from> FromType;
            typedef EntityRep<codim_to> ToType;
            typedef OrientedEntityRange<codim_to> row_type; // ??? doxygen henter doc fra Opm::SparseTable
            typedef Opm::SparseTable<ToType> super_t;
            typedef typename super_t::mutable_row_type mutable_row_type;

            OrientedEntityTable()
            {
            }

            template <typename DataIter, typename IntegerIter>
            OrientedEntityTable(DataIter data_beg, DataIter data_end,
                                IntegerIter rowsize_beg, IntegerIter rowsize_end)
                : super_t(data_beg, data_end, rowsize_beg, rowsize_end)
            {
            }

            using super_t::empty;
            using super_t::size;
            using super_t::dataSize;
            using super_t::clear;
            using super_t::appendRow;
            using super_t::allocate;

            int rowSize(const FromType& e) const
            {
                return super_t::rowSize(e.index());
            }

            row_type operator[](const FromType& e) const
            {
                return row_type(super_t::operator[](e.index()), e.orientation());
            }

            mutable_row_type row(const FromType& e)
            {
                return super_t::operator[](e.index());
            }
            bool operator==(const OrientedEntityTable& other) const
            {
                return super_t::operator==(other);
            }

            void swap(OrientedEntityTable& other)
            {
                super_t::swap(other);
            }

            void printSparseRelationMatrix(std::ostream& os) const
            {
                for (int i = 0; i < size(); ++i) {
                    const FromType from_ent(i, true);
                    const row_type r = operator[](from_ent);
                    const int rsize = r.size();
                    for (int j = 0; j < rsize; ++j) {
                        os << i << ' ' << r[j].index() << ' ' << (r[j].orientation() ? 1 : -1) << '\n';
                    }
                }
                os << std::flush;
            }

            void printRelationMatrix(std::ostream& os) const
            {
                int columns = numberOfColumns();
                for (int i = 0; i < size(); ++i) {
                    FromType from_ent(i, true);
                    row_type r  = operator[](from_ent);
                    int cur_col = 0;
                    int next_ent = 0;
                    ToType to_ent = r[next_ent];
                    int next_print = to_ent.index();
                    while (cur_col < columns) {
                        if (cur_col == next_print) {
                            if (to_ent.orientation()) {
                                os << "  1";
                            } else {
                                os << " -1";
                            }
                            ++next_ent;
                            if (next_ent >= r.size()) {
                                next_print = columns;
                            } else {
                                to_ent = r[next_ent];
                                next_print = to_ent.index();
                            }
                        } else {
                            os << "  0";
                        }
                        ++cur_col;
                    }
                    os << '\n';
                }
            }

            void makeInverseRelation(OrientedEntityTable<codim_to, codim_from>& inv) const
            {
                // Find the maximum index used. This will give (one less than) the size
                // of the table to be created.
                int maxind = -1;
                for (int i = 0; i < size(); ++i) {
                    EntityRep<codim_from> from_ent(i, true);
                    row_type r = operator[](from_ent);
                    for (int j = 0; j < r.size(); ++j) {
                        EntityRep<codim_to> to_ent = r[j];
                        int ind = to_ent.index();
                        maxind = std::max(ind, maxind);
                    }
                }
                // Build the new_sizes vector and compute datacount.
                std::vector<int> new_sizes(maxind + 1);
                int datacount = 0;
                for (int i = 0; i < size(); ++i) {
                    EntityRep<codim_from> from_ent(i, true);
                    row_type r = operator[](from_ent);
                    datacount += r.size();
                    for (int j = 0; j < r.size(); ++j) {
                        EntityRep<codim_to> to_ent = r[j];
                        int ind = to_ent.index();
                        ++new_sizes[ind];
                    }
                }
                // Compute the cumulative sizes.
                std::vector<int> cumul_sizes(new_sizes.size() + 1);
                cumul_sizes[0] = 0;
                std::partial_sum(new_sizes.begin(), new_sizes.end(), cumul_sizes.begin() + 1);
                // Using the cumulative sizes array as indices, we populate new_data.
                // Note that cumul_sizes[ind] is not kept constant, but incremented so that
                // it always gives the correct index for new data corresponding to index ind.
                std::vector<EntityRep<codim_from> > new_data(datacount);
                for (int i = 0; i < size(); ++i) {
                    EntityRep<codim_from> from_ent(i, true);
                    row_type r = operator[](from_ent);
                    for (int j = 0; j < r.size(); ++j) {
                        EntityRep<codim_to> to_ent(r[j]);
                        int ind = to_ent.index();
                        int data_ind = cumul_sizes[ind];
                        new_data[data_ind] = to_ent.orientation() ? from_ent : from_ent.opposite();
                        ++cumul_sizes[ind];
                    }
                }
                inv = OrientedEntityTable<codim_to, codim_from>(new_data.begin(),
                                                                new_data.end(),
                                                                new_sizes.begin(),
                                                                new_sizes.end());
            }

        private:
            int numberOfColumns() const
            {
                int maxind = 0;
                for (int i = 0; i < size(); ++i) {
                    FromType from_ent(i, true);
                    row_type r  = operator[](from_ent);
                    for (int j = 0; j < r.size(); ++j) {
                        maxind = std::max(maxind, r[j].index());
                    }
                }
                return maxind + 1;
            }
        };


    } // namespace cpgrid
} // namespace Dune




#endif // OPM_ORIENTEDENTITYTABLE_HEADER