linteractionvolume3d.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_FVMPFAL3DINTERACTIONVOLUME_HH
#define DUMUX_FVMPFAL3DINTERACTIONVOLUME_HH
 
#include "properties.hh"
 
namespace Dumux
{
 
// Mapper for local interaction volume indices (see doc/docextra/3dmpfa).
class IndexTranslator
{
public:
    enum
    {
        subVolumeTotalNum = 8,
        fluxFacesTotalNum = 12,
        fluxFacesNumOnSubVolume = 3,
        fluxEdgesTotalNum = 6,
        edgesNumOnFluxFace = 2
    };
 
    static int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
    }
 
    static int getEdgeIndexFromSubVolumeFace(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
    {
        return edgeIndexOnSubVolumeFace_[subVolumeIdx][subVolumeFaceIdx][subVolumeEdgeIdx];
    }
 
    static int getFaceIndexFromElements(int elem1Idx, int elem2Idx)
    {
        return faceIndexFromElements_[elem1Idx][elem2Idx];
    }
 
 
private:
    static const int faceIndexOnSubVolume_[subVolumeTotalNum][fluxFacesNumOnSubVolume];
    static const int edgeIndexOnSubVolumeFace_[subVolumeTotalNum][fluxFacesNumOnSubVolume][edgesNumOnFluxFace];
    static const int faceIndexFromElements_[subVolumeTotalNum][subVolumeTotalNum];
};
 
const int IndexTranslator::faceIndexOnSubVolume_[subVolumeTotalNum][fluxFacesNumOnSubVolume] =
{
        {0, 3, 8},
        {1, 0, 9},
        {3, 2, 11},
        {2, 1, 10},
        {8, 4, 7},
        {9, 5, 4},
        {11, 7, 6},
        {10, 6, 5}
};
 
const int IndexTranslator::edgeIndexOnSubVolumeFace_[subVolumeTotalNum][fluxFacesNumOnSubVolume][edgesNumOnFluxFace] =
{
        {{2, 0},{0, 5},{5, 2}},
        {{3, 0},{0, 2},{2, 3}},
        {{5, 0},{0, 4},{4, 5}},
        {{4, 0},{0, 3},{3, 4}},
        {{2, 5},{1, 2},{5, 1}},
        {{3, 2},{1, 3},{2, 1}},
        {{5, 4},{1, 5},{4, 1}},
        {{4, 3},{1, 4},{3, 1}}
};
 
const int IndexTranslator::faceIndexFromElements_[subVolumeTotalNum][subVolumeTotalNum] =
{
        {-1, 0, 3, -1, 8, -1, -1, -1},
        {0, -1, -1, 1, -1, 9, -1, -1},
        {3, -1, -1, 2, -1, -1, 11, -1},
        {-1, 1, 2, -1, -1, -1, -1, 10},
        {8, -1, -1, -1, -1, 4, 7, -1},
        {-1, 9, -1, -1, 4, -1, -1, 5},
        {-1, -1, 11, -1, 7, -1, -1, 6},
        {-1, -1, -1, 10, -1, 5, 6, -1}
};
 
 
template<class TypeTag>
class FvMpfaL3dInteractionVolume
{
private:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
 
    enum
    {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,
    };
 
    using Element = typename GridView::template Codim<0>::Entity;
    using ElementSeed = typename Grid::template Codim<0>::EntitySeed;
 
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
    using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
    using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
    using DimVector = Dune::FieldVector<Scalar, dim>;
    using FieldVectorVector = Dune::FieldVector<DimVector, dim>;
    using FieldVectorVector2 = Dune::FieldVector<DimVector, 2>;
    using FieldVectorVectorVector = Dune::FieldVector<FieldVectorVector2, dim>;
    using IndexVector = Dune::FieldVector<int, dim>;
    using BCTypeVector = std::vector<BoundaryTypes>;
    using BCVector = std::vector<PrimaryVariables>;
 
public:
    //An interaction volume around a vertex includes in general 12 flux faces (see doc/docextra/3dmpfa).
    //If the vertex is on the boundary these flux faces can be inside the domain on the boundary or outside the domain.
    enum FaceTypes
    {
        inside = 1,
        boundary = 0,
        outside = -1,
    };
 
    enum
    {
        subVolumeTotalNum = IndexTranslator::subVolumeTotalNum,
        fluxFacesTotalNum = IndexTranslator::fluxFacesTotalNum,
        fluxEdgesTotalNum = IndexTranslator::fluxEdgesTotalNum
    };
 
 
    FvMpfaL3dInteractionVolume(const Grid& grid)
    : grid_(&grid)
    , normal_(FieldVectorVector(DimVector(0.0)))
    , facePos_(DimVector(0.0))
    , edgePos_((DimVector(0.0)))
    , faceArea_(0.0)
    , faceType_(fluxFacesTotalNum, inside)
    , indexOnElement_(IndexVector(0.0))
    , elements_(subVolumeTotalNum)
    , centerVertexPos_(0)
    , elementNum_(0)
    {}
 
    void reset()
    {
        elements_.clear();
        elements_.resize(subVolumeTotalNum);
        centerVertexPos_ = 0;
        indexOnElement_ = IndexVector(0.0);
        faceType_.clear();
        faceType_.resize(fluxFacesTotalNum, inside);
        faceArea_ = 0.0;
        facePos_ = DimVector(0.0);
        edgePos_ = DimVector(0.0);
        normal_ = FieldVectorVector(DimVector(0.0));
        elementNum_ = 0;
        boundaryTypes_.clear();
        neumannValues_.clear();
        dirichletValues_.clear();
    }
 
    /*
     * \param centerVertexPos Position of the central vertex
     */
    void setCenterPosition(const DimVector &centerVertexPos)
    {
        centerVertexPos_ = centerVertexPos;
    }
 
    void setSubVolumeElement(const Element& element, int subVolumeIdx)
    {
        if (!hasSubVolumeElement(subVolumeIdx))
        {
            elements_[subVolumeIdx].push_back(element.seed());
            elementNum_++;
        }
        else
        {
            elements_[subVolumeIdx].insert(elements_[subVolumeIdx].begin(), element.seed());
        }
    }
 
 
    void setFacePosition(const DimVector& pos, int fIdx)
    {
        facePos_[fIdx] = pos;
    }
 
 
    void setEdgePosition(const DimVector& pos, int edgeIdx)
    {
        edgePos_[edgeIdx] = pos;
    }
 
 
    void setFaceArea(Scalar faceArea, int fIdx)
    {
        faceArea_[fIdx] = faceArea;
    }
 
 
    void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdx)
    {
        normal_[subVolumeIdx][subVolumeFaceIdx] = normal;
    }
 
 
    void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
    {
        if (boundaryTypes_.size() == 0)
        {
            boundaryTypes_.resize(fluxFacesTotalNum);
        }
        boundaryTypes_[subVolumeFaceIdx] = boundaryTypes;
        faceType_[subVolumeFaceIdx] = boundary;
    }
 
 
    void setOutsideFace(int subVolumeFaceIdx)
    {
        faceType_[subVolumeFaceIdx] = outside;
    }
 
 
    void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
    {
        if (dirichletValues_.size() == 0)
        {
            dirichletValues_.resize(fluxFacesTotalNum, PrimaryVariables(0.0));
        }
        dirichletValues_[subVolumeFaceIdx] = condition;
    }
 
 
    void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
    {
        if (neumannValues_.size() == 0)
        {
            neumannValues_.resize(fluxFacesTotalNum, PrimaryVariables(0.0));
        }
        neumannValues_[subVolumeFaceIdx] = condition;
    }
 
 
    void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdx)
    {
        indexOnElement_[subVolumeIdx][subVolumeFaceIdx] = indexInInside;
    }
 
    DimVector& getCenterPosition()
    {
        return centerVertexPos_;
    }
 
    /*
     * \param subVolumeIdx The local element index in the interaction volume
     * \param subVolumeFaceIdx The local face index in the interaction volume element
     *
     * \return The interaction volume face index
     */
    int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx);
    }
 
    /*
     * \param subVolumeIdx The local element index in the interaction volume
     * \param subVolumeFaceIdx The local face index in the interaction volume element
     * \param subVolumeEdgeIdx The local edge index in the interaction volume element
     *
     * \return The interaction volume edge index
     */
    int getEdgeIndexFromSubVolumeFace(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
    {
        return IndexTranslator::getEdgeIndexFromSubVolumeFace(subVolumeIdx, subVolumeFaceIdx, subVolumeEdgeIdx);
    }
 
    int getElementNumber()
    {
        return elementNum_;
    }
 
 
    int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
    }
 
 
    Element getSubVolumeElement(int subVolumeIdx)
    {
        if (hasSubVolumeElement(subVolumeIdx))
            return grid_->entity(elements_[subVolumeIdx][0]);
        else
        {
            std::cout<<"Problems when calling getSubVolumeElement("<<subVolumeIdx<<")\n";
            printInteractionVolumeInfo();
            DUNE_THROW(Dune::RangeError, "element not in interaction volume!");
        }
    }
 
 
    bool hasSubVolumeElement(int subVolumeIdx)
    {
        return (elements_[subVolumeIdx].size() > 0);
    }
 
 
    BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
    {
        return boundaryTypes_[subVolumeFaceIdx];
    }
 
 
    bool isOutsideFace(int subVolumeFaceIdx)
    {
        return (faceType_[subVolumeFaceIdx] == outside);
    }
 
 
    bool isInsideFace(int subVolumeFaceIdx)
    {
        return (faceType_[subVolumeFaceIdx] == inside);
    }
 
 
    bool isBoundaryFace(int subVolumeFaceIdx)
    {
        return (faceType_[subVolumeFaceIdx] == boundary);
    }
 
 
    bool isInnerVolume()
    {
        for (unsigned int i = 0; i < faceType_.size(); i++)
        {
            if (isOutsideFace(i) || isBoundaryFace(i))
                return false;
        }
        return true;
    }
 
 
    bool isBoundaryInteractionVolume()
    {
        return (boundaryTypes_.size() > 0);
    }
 
 
    PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
    {
        return dirichletValues_[subVolumeFaceIdx];
    }
 
 
    PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
    {
        return neumannValues_[subVolumeFaceIdx];
    }
 
 
    DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return normal_[subVolumeIdx][subVolumeFaceIdx];
    }
 
 
    DimVector& getFacePosition(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return facePos_[IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx)];
    }
 
 
    DimVector& getEdgePosition(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
    {
        return edgePos_[IndexTranslator::getEdgeIndexFromSubVolumeFace(subVolumeIdx, subVolumeFaceIdx, subVolumeEdgeIdx)];
    }
 
 
    Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdx)
    {
        return faceArea_[IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx)];
    }
 
 
    DimVector& getFacePosition(int fIdx)
    {
        return facePos_[fIdx];
    }
 
 
    DimVector& getEdgePosition(int edgeIdx)
    {
        return edgePos_[edgeIdx];
    }
 
 
    Scalar& getFaceArea(int fIdx)
    {
        return faceArea_[fIdx];
    }
 
    void printInteractionVolumeInfo()
    {
        std::cout<<"\nNumber of stored elements: "<<elementNum_<<"\n";
        std::cout<<" center position: "<<centerVertexPos_<<"\n";
        for (int i = 0; i < subVolumeTotalNum; i++)
        {
            if (elements_[i].size() > 0)
            {
            std::cout<<"element "<<i<<":\n";
            std::cout<<"element level: "<<grid_->entity(elements_[i][0]).level()<<"\n";
            std::cout<<"element position: "<<grid_->entity(elements_[i][0]).geometry().center()<<"\n";
            std::cout<<"element volume: "<<grid_->entity(elements_[i][0]).geometry().volume()<<"\n";
            std::cout<<"face indices on element: "<<indexOnElement_[i]<<"\n";
            std::cout<<"face normals on element: "<<normal_[i]<<"\n";
            std::cout<<"face areas on element: ";
            for (int j = 0; j < 3; j++)
            {
                std::cout<<getFaceArea(i, j)<<" ";
            }
            std::cout<<"\n";
            std::cout<<"face position on element: ";
            for (int j = 0; j < 3; j++)
            {
                std::cout<<getFacePosition(i, j)<<" ";
            }
            std::cout<<"\n";
            std::cout<<"edgePos:";
            for (int j = 0; j < 3; j++)
            {
                std::cout<<"    ";
                for (int k = 0; k < 2; k++)
                {
                    std::cout<<getEdgePosition(i, j, k) <<" ";
                }
            }
            std::cout<<"\n";
            }
        }
    }
 
//    void printInteractionVolumeInfoToFile(std::ofstream& dataFile)
//    {
//        dataFile<<"\nNumber of stored elements: "<<elementNum_<<"\n";
//        dataFile<<" center position: "<<centerVertexPos_<<"\n";
//        for (int i = 0; i < subVolumeTotalNum; i++)
//        {
//            if (elements_[i].size() > 0)
//            {
//            dataFile<<"element "<<i<<":\n";
//            dataFile<<"element level: "<<elements_[i][0]->level()<<"\n";
//            dataFile<<"element position: "<<elements_[i][0]->geometry().center()<<"\n";
//            dataFile<<"element volume: "<<elements_[i][0]->geometry().volume()<<"\n";
//            dataFile<<"face indices on element: "<<indexOnElement_[i]<<"\n";
//            dataFile<<"face normals on element: "<<normal_[i]<<"\n";
//            dataFile<<"face areas on element: ";
//            for (int j = 0; j < 3; j++)
//            {
//                dataFile<<getFaceArea(i, j)<<" ";
//            }
//            dataFile<<"\n";
//            dataFile<<"face position on element: ";
//            for (int j = 0; j < 3; j++)
//            {
//                dataFile<<getFacePosition(i, j)<<" ";
//            }
//            dataFile<<"\n";
//            dataFile<<"edgePos:";
//            for (int j = 0; j < 3; j++)
//            {
//                dataFile<<"    ";
//                for (int k = 0; k < 2; k++)
//                {
//                    dataFile<<getEdgePosition(i, j, k) <<" ";
//                }
//            }
//            dataFile<<"\n";
//
//        }
//        dataFile<<"face types: ";
//        for (int i = 0; i < fluxFacesTotalNum; i++)
//        {
//            dataFile<<faceType_[i]<<" ";
//        }
//        dataFile<<"\n";
//        if (isBoundaryInteractionVolume())
//        {
//        dataFile<<"Boundaries:\n";
//        dataFile<<"dirichlet: ";
//        for (int i = 0; i < fluxFacesTotalNum; i++)
//        {
//            dataFile<<boundaryTypes_[i].isDirichlet(0)<<" ";
//            dataFile<<boundaryTypes_[i].isDirichlet(1)<<" ";
//        }
//        dataFile<<"\n";
//        dataFile<<"neumann: ";
//        for (int i = 0; i < fluxFacesTotalNum; i++)
//        {
//            dataFile<<boundaryTypes_[i].isNeumann(0)<<" ";
//            dataFile<<boundaryTypes_[i].isNeumann(1)<<" ";
//        }
//        dataFile<<"\n";
//        dataFile<<"values: ";
//        for (int i = 0; i < fluxFacesTotalNum; i++)
//        {
//            if (dirichletValues_.size() > 0)
//            dataFile<<dirichletValues_[i]<<" ";
//            if (neumannValues_.size() > 0)
//            dataFile<<neumannValues_[i]<<" ";
//        }
//        dataFile<<"\n";
//        }
//    }
 
private:
    const Grid* grid_;
    Dune::FieldVector<FieldVectorVector, subVolumeTotalNum> normal_;
    Dune::FieldVector<DimVector, fluxFacesTotalNum> facePos_;
    Dune::FieldVector<DimVector, fluxEdgesTotalNum> edgePos_;
    Dune::FieldVector<Scalar, fluxFacesTotalNum> faceArea_;
    BCTypeVector boundaryTypes_;
    std::vector<int> faceType_;
    Dune::FieldVector<IndexVector, subVolumeTotalNum> indexOnElement_;
    std::vector<std::vector<ElementSeed> > elements_;
    BCVector neumannValues_;
    BCVector dirichletValues_;
    DimVector centerVertexPos_;
    int elementNum_;
};
}
#endif