facecentered/staggered/connectivitymap.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:
/*****************************************************************************
 *   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/>.   *
 *****************************************************************************/
#ifndef DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_CONNECTIVITY_MAP_HH
#define DUMUX_DISCRETIZATION_FACECENTERED_STAGGERED_CONNECTIVITY_MAP_HH
 
#include <algorithm>
#include <vector>
#include <dune/grid/common/partitionset.hh>
#include <dune/grid/common/gridenums.hh>
#include <dumux/common/indextraits.hh>
 
namespace Dumux {
 
template<class GridGeometry>
class FaceCenteredStaggeredConnectivityMap
{
    using GridView = typename GridGeometry::GridView;
    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
    using Stencil = std::vector<GridIndexType>;
    using Map = std::vector<Stencil>;
 
public:
 
    void update(const GridGeometry& gridGeometry)
    {
        map_.clear();
        map_.resize(gridGeometry.numScv());
 
        auto fvGeometry = localView(gridGeometry);
        for (const auto& element : elements(gridGeometry.gridView()))
        {
            if (element.partitionType() == Dune::InteriorEntity)
                continue;
 
            assert(element.partitionType() == Dune::OverlapEntity);
 
            // restrict the FvGeometry locally and bind to the element
            fvGeometry.bind(element);
            // loop over sub control faces
            for (const auto& scvf : scvfs(fvGeometry))
            {
                if (scvf.isFrontal() && !scvf.boundary() && !scvf.processorBoundary())
                {
                    const auto& ownScv = fvGeometry.scv(scvf.insideScvIdx());
                    const auto& facet = element.template subEntity <1> (ownScv.indexInElement());
                    if (facet.partitionType() == Dune::BorderEntity)
                        map_[ownScv.index()].push_back(scvf.outsideScvIdx());
                }
            }
        }
 
        for (const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
        {
            fvGeometry.bind(element);
 
            // loop over sub control faces
            for (const auto& scvf : scvfs(fvGeometry))
            {
                assert(!scvf.processorBoundary());
                const auto& ownScv = fvGeometry.scv(scvf.insideScvIdx());
                const auto ownDofIndex = ownScv.dofIndex();
                const auto ownScvIndex = ownScv.index();
 
                if (scvf.isFrontal())
                {
                    if (!scvf.boundary()) // opposite dof
                        map_[ownScvIndex].push_back(scvf.outsideScvIdx());
                    else
                    {
                        // treat frontal faces on boundaries
                        for (const auto& scv : scvs(fvGeometry))
                        {
                            const auto otherDofIndex = scv.dofIndex();
                            if (ownDofIndex != otherDofIndex)
                                map_[scv.index()].push_back(ownScvIndex);
                        }
                    }
                }
                else // lateral face
                {
                    // the parallel DOF scv
                    // outsideScv   insideScv
                    // ###y######|s|#####x####
                    //           |c|
                    //           |v|
                    //           |f|
                    if (!scvf.boundary())
                        map_[ownScvIndex].push_back(scvf.outsideScvIdx());
 
 
                    // the normal DOF scv
                    //
                    //     outsideScv
                    //
                    // |s|orthogonalScvf
                    // |c|
                    // |v| insideScv
                    // |f|
                    const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
                    assert(orthogonalScvf.isLateral());
                    map_[ownScvIndex].push_back(orthogonalScvf.insideScvIdx());
                    if (!orthogonalScvf.boundary())
                        map_[ownScvIndex].push_back(orthogonalScvf.outsideScvIdx());
                }
            }
        }
 
        // make stencils unique
        for (auto& stencil : map_)
        {
            std::sort(stencil.begin(), stencil.end());
            stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
        }
    }
 
    const Stencil& operator[] (const GridIndexType globalI) const
    { return map_[globalI]; }
 
private:
    Map map_;
};
 
} // end namespace Dumux
 
#endif