Commit 692aee82 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[porenetworkmodel] Add files for 1p, 1pnc and 2p PNM

parent 1f613dab
......@@ -12,6 +12,7 @@ add_subdirectory(material)
add_subdirectory(multidomain)
add_subdirectory(nonlinear)
add_subdirectory(parallel)
add_subdirectory(porenetwork)
add_subdirectory(porousmediumflow)
# if Python bindings are enabled, include necessary sub directories.
......
add_subdirectory(box)
add_subdirectory(cellcentered)
add_subdirectory(fem)
add_subdirectory(porenetwork)
add_subdirectory(projection)
add_subdirectory(staggered)
......
install(FILES
fvelementgeometry.hh
gridgeometry.hh
subcontrolvolume.hh
subcontrolvolumeface.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/porenetwork)
// -*- 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/>. *
*****************************************************************************/
/*!
* \file
* \ingroup PoreNetworkModels
* \brief Base class for the local geometry for porenetworks
*/
#ifndef DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/scvandscvfiterators.hh>
namespace Dumux::PoreNetwork {
/*!
* \ingroup PoreNetworkModels
* \brief Base class for the local geometry for porenetworks
* \tparam GG the finite volume grid geometry type
* \tparam enableFVGridGeometryCache if the grid geometry is cached or not
*/
template<class GG, bool enableFVGridGeometryCache>
class PNMFVElementGeometry;
//! specialization in case the FVElementGeometries are stored
template<class GG>
class PNMFVElementGeometry<GG, true>
{
using GridView = typename GG::GridView;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using Element = typename GridView::template Codim<0>::Entity;
using CoordScalar = typename GridView::ctype;
using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
public:
//! export type of subcontrol volume
using SubControlVolume = typename GG::SubControlVolume;
//! export type of subcontrol volume face
using SubControlVolumeFace = typename GG::SubControlVolumeFace;
//! export type of finite volume grid geometry
using GridGeometry = GG;
//! the maximum number of scvs per element
static constexpr std::size_t maxNumElementScvs = 2;
//! Constructor
PNMFVElementGeometry(const GridGeometry& gridGeometry)
: gridGeometryPtr_(&gridGeometry) {}
//! Get a sub control volume with a local scv index
const SubControlVolume& scv(LocalIndexType scvIdx) const
{
return gridGeometry().scvs(eIdx_)[scvIdx];
}
//! Get a sub control volume face with a local scvf index
const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
{
return gridGeometry().scvfs(eIdx_)[scvfIdx];
}
//! iterator range for sub control volumes. Iterates over
//! all scvs of the bound element.
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
{
const auto& g = fvGeometry.gridGeometry();
using Iter = typename std::decay_t<decltype(g.scvs(fvGeometry.eIdx_))>::const_iterator;
return Dune::IteratorRange<Iter>(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end());
}
//! iterator range for sub control volumes faces. Iterates over
//! all scvfs of the bound element.
//! This is a free function found by means of ADL
//! To iterate over all sub control volume faces of this FVElementGeometry use
//! for (auto&& scvf : scvfs(fvGeometry))
friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
{
const auto& g = fvGeometry.gridGeometry();
using Iter = typename std::decay_t<decltype(g.scvfs(fvGeometry.eIdx_))>::const_iterator;
return Dune::IteratorRange<Iter>(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end());
}
//! Get a local finite element basis
const FeLocalBasis& feLocalBasis() const
{
return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
}
//! The total number of sub control volumes
std::size_t numScv() const
{
return gridGeometry().scvs(eIdx_).size();
}
//! The total number of sub control volume faces
std::size_t numScvf() const
{
return gridGeometry().scvfs(eIdx_).size();
}
//! this function is for compatibility reasons with cc methods
//! The box stencil is always element-local so bind and bindElement
//! are identical.
void bind(const Element& element)
{
this->bindElement(element);
}
//! Binding of an element, has to be called before using the fvgeometries
//! Prepares all the volume variables within the element
//! For compatibility reasons with the FVGeometry cache being disabled
void bindElement(const Element& element)
{
elementPtr_ = &element;
eIdx_ = gridGeometry().elementMapper().index(element);
}
//! The global finite volume geometry we are a restriction of
const GridGeometry& gridGeometry() const
{ return *gridGeometryPtr_; }
//! Returns whether one of the geometry's scvfs lies on a boundary
bool hasBoundaryScvf() const
{ return gridGeometry().hasBoundaryScvf(eIdx_); }
private:
const Element* elementPtr_;
const GridGeometry* gridGeometryPtr_;
GridIndexType eIdx_;
};
//! specialization in case the FVElementGeometries are not stored
template<class GG>
class PNMFVElementGeometry<GG, false>
{
using GridView = typename GG::GridView;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using Element = typename GridView::template Codim<0>::Entity;
using CoordScalar = typename GridView::ctype;
using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
public:
//! export type of subcontrol volume
using SubControlVolume = typename GG::SubControlVolume;
//! export type of subcontrol volume face
using SubControlVolumeFace = typename GG::SubControlVolumeFace;
//! export type of finite volume grid geometry
using GridGeometry = GG;
//! the maximum number of scvs per element
static constexpr std::size_t maxNumElementScvs = 2;
//! Constructor
PNMFVElementGeometry(const GridGeometry& gridGeometry)
: gridGeometryPtr_(&gridGeometry) {}
//! Get a sub control volume with a local scv index
const SubControlVolume& scv(LocalIndexType scvIdx) const
{
return scvs_[scvIdx];
}
//! Get a sub control volume face with a local scvf index
const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
{
return scvfs_[0];
}
//! iterator range for sub control volumes. Iterates over
//! all scvs of the bound element.
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
{
using Iter = typename std::decay_t<decltype(fvGeometry.scvs_)>::const_iterator;
return Dune::IteratorRange<Iter>(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end());
}
//! iterator range for sub control volumes faces. Iterates over
//! all scvfs of the bound element.
//! This is a free function found by means of ADL
//! To iterate over all sub control volume faces of this FVElementGeometry use
//! for (auto&& scvf : scvfs(fvGeometry))
friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
{
using Iter = typename std::decay_t<decltype(fvGeometry.scvfs_)>::const_iterator;
return Dune::IteratorRange<Iter>(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end());
}
//! Get a local finite element basis
const FeLocalBasis& feLocalBasis() const
{
return gridGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
}
//! The total number of sub control volumes
std::size_t numScv() const
{
return scvs_.size();
}
//! The total number of sub control volume faces
std::size_t numScvf() const
{
return scvfs_.size();
}
//! this function is for compatibility reasons with cc methods
//! The box stencil is always element-local so bind and bindElement
//! are identical.
void bind(const Element& element)
{
this->bindElement(element);
}
//! Binding of an element, has to be called before using the fvgeometries
//! Prepares all the volume variables within the element
//! For compatibility reasons with the FVGeometry cache being disabled
void bindElement(const Element& element)
{
elementPtr_ = &element;
eIdx_ = gridGeometry().elementMapper().index(element);
makeElementGeometries(element);
}
//! The global finite volume geometry we are a restriction of
const GridGeometry& gridGeometry() const
{ return *gridGeometryPtr_; }
//! Returns whether one of the geometry's scvfs lies on a boundary
bool hasBoundaryScvf() const
{ return hasBoundaryScvf_; }
private:
void makeElementGeometries(const Element& element)
{
hasBoundaryScvf_ = false;
// get the element geometry
auto elementGeometry = element.geometry();
// construct the sub control volumes
for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
{
// get asssociated dof index
const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
// get the corners
auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()};
// get the fractional volume asssociated with the scv
const auto volume = gridGeometry().poreVolume(dofIdxGlobal) / gridGeometry().coordinationNumber(dofIdxGlobal);
// add scv to the local container
scvs_[scvLocalIdx] = SubControlVolume(dofIdxGlobal,
scvLocalIdx,
eIdx_,
std::move(corners),
volume);
if (gridGeometry().poreLabel(dofIdxGlobal) > 0)
hasBoundaryScvf_ = true;
}
// construct the inner sub control volume face
auto unitOuterNormal = elementGeometry.corner(1) - elementGeometry.corner(0);
unitOuterNormal /= unitOuterNormal.two_norm();
LocalIndexType scvfLocalIdx = 0;
scvfs_[0] = SubControlVolumeFace(elementGeometry.center(),
std::move(unitOuterNormal),
gridGeometry().throatCrossSectionalArea(gridGeometry().elementMapper().index(element)),
scvfLocalIdx++,
std::vector<LocalIndexType>({0, 1}));
}
//! The bound element
const Element* elementPtr_;
GridIndexType eIdx_;
//! The global geometry this is a restriction of
const GridGeometry* gridGeometryPtr_;
//! vectors to store the geometries locally after binding an element
std::array<SubControlVolume, 2> scvs_;
std::array<SubControlVolumeFace, 1> scvfs_;
bool hasBoundaryScvf_ = false;
};
} // end namespace Dumux::PoreNetwork
#endif
This diff is collapsed.
// -*- 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/>. *
*****************************************************************************/
/*!
* \file
* \ingroup PoreNetworkModels
* \brief the sub control volume for pore networks
*/
#ifndef DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUME_HH
#define DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUME_HH
#include <dune/geometry/affinegeometry.hh>
#include <dumux/common/math.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/subcontrolvolumebase.hh>
namespace Dumux::PoreNetwork {
/*!
* \ingroup PoreNetworkModels
* \brief Default traits class
* \tparam GV the type of the grid view
*/
template<class GridView>
struct PNMDefaultScvGeometryTraits
{
using Grid = typename GridView::Grid;
static const int dim = Grid::dimension;
static const int dimWorld = Grid::dimensionworld;
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using Scalar = typename Grid::ctype;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CornerStorage = std::array<GlobalPosition, 2>;
using Geometry = Dune::AffineGeometry<Scalar, 1, dimWorld>;
};
/*!
* \ingroup PoreNetworkModels
* \brief the sub control volume for porenetworks
* \tparam GV the type of the grid view
* \tparam T the scv geometry traits
*/
template<class GV,
class T = PNMDefaultScvGeometryTraits<GV> >
class PNMSubControlVolume
: public Dumux::SubControlVolumeBase<PNMSubControlVolume<GV, T>, T>
{
using ThisType = PNMSubControlVolume<GV, T>;
using ParentType = Dumux::SubControlVolumeBase<ThisType, T>;
using GridIndexType = typename T::GridIndexType;
using LocalIndexType = typename T::LocalIndexType;
using Scalar = typename T::Scalar;
using CornerStorage = typename T::CornerStorage;
using Geometry = typename T::Geometry;
public:
//! export the type used for global coordinates
using GlobalPosition = typename T::GlobalPosition;
//! state the traits public and thus export all types
using Traits = T;
//! The default constructor
PNMSubControlVolume() = default;
// the contructor in the box case
template<class Corners>
PNMSubControlVolume(GridIndexType dofIndex,
LocalIndexType scvIdx,
GridIndexType elementIndex,
Corners&& corners,
const Scalar volume)
: center_((corners[0]+corners[1])/2.0),
corners_(std::forward<CornerStorage>(corners)),
volume_(volume),
elementIndex_(elementIndex),
localDofIdx_(scvIdx),
dofIndex_(dofIndex)
{}
//! The center of the sub control volume (return pore center).
//! Be aware that this is not the pore-body center! Use dofPosition() for the latter!
const GlobalPosition& center() const
{ return center_; }
//! The volume of the sub control volume (part of a pore)
Scalar volume() const
{ return volume_; }
//! The element-local index of the dof this scv is embedded in
LocalIndexType localDofIndex() const
{ return localDofIdx_; }
//! The element-local index of this scv.
LocalIndexType indexInElement() const
{ return localDofIdx_; }
//! The index of the dof this scv is embedded in
GridIndexType dofIndex() const
{ return dofIndex_; }
// The position of the dof this scv is embedded in
const GlobalPosition& dofPosition() const
{ return corners_[0]; }
//! The global index of the element this scv is embedded in
GridIndexType elementIndex() const
{ return elementIndex_; }
//! Return the corner for the given local index
const GlobalPosition& corner(LocalIndexType localIdx) const
{
assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
return corners_[localIdx];
}
//! The geometry of the sub control volume e.g. for integration
Geometry geometry() const
{
return Geometry(Dune::GeometryTypes::simplex(1), corners_);
}
private:
GlobalPosition center_;
CornerStorage corners_;
Scalar volume_;
GridIndexType elementIndex_;
LocalIndexType localDofIdx_;
GridIndexType dofIndex_;
};
} // end namespace Dumux::PoreNetwork
#endif
// -*- 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/>. *
*****************************************************************************/
/*!
* \file
* \ingroup PoreNetworkModels
* \brief Base class for a sub control volume face
*/
#ifndef DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUMEFACE_HH
#define DUMUX_DISCRETIZATION_PNM_SUBCONTROLVOLUMEFACE_HH
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/subcontrolvolumefacebase.hh>
namespace Dumux::PoreNetwork {
/*!
* \ingroup PoreNetworkModels
* \brief Default traits class
* \tparam GV the type of the grid view
*/
template<class GridView>
struct PNMDefaultScvfGeometryTraits
{
using Grid = typename GridView::Grid;
static constexpr int dim = Grid::dimension;
static constexpr int dimWorld = Grid::dimensionworld;
using GridIndexType = typename IndexTraits<GridView>::GridIndex;
using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
using Scalar = typename Grid::ctype;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CornerStorage = std::array<Dune::FieldVector<Scalar, dimWorld>, 1>;
};
/*!
* \ingroup PoreNetworkDiscretization
* \brief Class for a sub control volume face for porenetworks
* \tparam GV the type of the grid view
* \tparam T the scvf geometry traits
*/
template<class GV,
class T = PNMDefaultScvfGeometryTraits<GV> >
class PNMSubControlVolumeFace
: public Dumux::SubControlVolumeFaceBase<PNMSubControlVolumeFace<GV, T>, T>
{
using ThisType = PNMSubControlVolumeFace<GV, T>;
using ParentType = Dumux::SubControlVolumeFaceBase<ThisType, T>;
using GridIndexType = typename T::GridIndexType;
using LocalIndexType = typename T::LocalIndexType;
using Scalar = typename T::Scalar;
public:
//! export the type used for global coordinates
using GlobalPosition = typename T::GlobalPosition;
//! state the traits public and thus export all types
using Traits = T;
//! The default constructor
PNMSubControlVolumeFace() = default;
//! Constructor for inner scvfs
PNMSubControlVolumeFace(const GlobalPosition& center,
const GlobalPosition& unitOuterNormal,
const Scalar& area,
GridIndexType scvfIndex,
std::vector<LocalIndexType>&& scvIndices)
: center_(center),
unitOuterNormal_(unitOuterNormal),
area_(area),
scvfIndex_(scvfIndex),
scvIndices_({scvIndices[0], scvIndices[1]})
{}
//! The center of the sub control volume face
const GlobalPosition& center() const
{ return center_; }
//! The integration point for flux evaluations in global coordinates
const GlobalPosition& ipGlobal() const