Commit adf2f340 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

add box-dfm framework

parent 53aec264
......@@ -109,7 +109,7 @@ public:
return 1.0;
}
private:
protected:
const typename Element::Geometry& elementGeometry_; //!< Reference to the element geometry
std::size_t corners_; // number of element corners
std::array<GlobalPosition, maxPoints> p; // the points needed for construction of the geometries
......@@ -337,7 +337,7 @@ public:
return (p[1]-p[0]).two_norm();
}
private:
protected:
const typename Element::Geometry& elementGeometry_; //!< Reference to the element geometry
std::size_t corners_; // number of element corners
std::array<GlobalPosition, maxPoints> p; // the points needed for construction of the geometries
......@@ -624,7 +624,7 @@ public:
return 0.5*Dumux::crossProduct(p[3]-p[0], p[2]-p[1]).two_norm();
}
private:
protected:
const typename Element::Geometry& elementGeometry_; //!< Reference to the element geometry
std::size_t corners_; // number of element corners
std::array<GlobalPosition, maxPoints> p; // the points needed for construction of the scv/scvf geometries
......
......@@ -9,6 +9,7 @@ add_subdirectory("2pncmin")
add_subdirectory("3p")
add_subdirectory("3pwateroil")
add_subdirectory("3p3c")
add_subdirectory("boxdfm")
add_subdirectory("co2")
add_subdirectory("compositional")
add_subdirectory("immiscible")
......
#install headers
install(FILES
fluxvariablescache.hh
fvelementgeometry.hh
fvgridgeometry.hh
geometryhelper.hh
model.hh
subcontrolvolume.hh
subcontrolvolumeface.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/boxdfm)
// -*- 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 2 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
* \brief Cache class for the flux variables to be used
* in conjunction with the box discrete fracture scheme
*/
#ifndef DUMUX_POROUSMEDIUM_BOXDFM_FLUXVARIABLESCACHE_HH
#define DUMUX_POROUSMEDIUM_BOXDFM_FLUXVARIABLESCACHE_HH
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/fluxvariablescaching.hh>
namespace Dumux {
/*!
* \ingroup BoxDiscretization
* \ingroup BoxDFM
* \brief We only store discretization-related quantities for the box method.
* However, we cannot reuse the cache of the standard box method as we have
* to take into account the scvs that lie on fracture facets.
*/
template<class TypeTag>
class BoxDfmFluxVariablesCache
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
using TransmissibilityVector = std::vector<IndexType>;
using CoordScalar = typename GridView::ctype;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
void update(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
{
const auto geometry = element.geometry();
const auto& localBasis = fvGeometry.feLocalBasis();
// evaluate shape functions and gradients at the integration point
std::vector<ShapeValue> shapeVals;
const auto ipLocal = geometry.local(scvf.ipGlobal());
jacInvT_ = geometry.jacobianInverseTransposed(ipLocal);
localBasis.evaluateJacobian(ipLocal, shapeJacobian_);
localBasis.evaluateFunction(ipLocal, shapeVals);
// set the shape values
shapeValues_.resize(fvGeometry.numScv(), 0.0);
if (!scvf.isOnFracture())
std::copy(shapeVals.begin(), shapeVals.end(), shapeValues_.begin());
else
{
const auto thisFacetIdx = scvf.facetIndexInElement();
for (const auto& scv: scvs(fvGeometry))
if (scv.isOnFracture() && scv.facetIndexInElement() == thisFacetIdx)
shapeValues_[scv.indexInElement()] = shapeVals[scv.localDofIndex()];
}
// set the shape value gradients
gradN_.resize(fvGeometry.numScv(), GlobalPosition(0.0));
if (!scvf.isOnFracture())
{
for (const auto& scv: scvs(fvGeometry))
if (!scv.isOnFracture())
jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
}
else
{
const auto thisFacetIdx = scvf.facetIndexInElement();
// first, find all local dofs on this facet
std::vector<unsigned int> facetLocalDofs;
for (const auto& scv : scvs(fvGeometry))
if (scv.isOnFracture() && scv.facetIndexInElement() == thisFacetIdx)
facetLocalDofs.push_back(scv.localDofIndex());
for (const auto& scv: scvs(fvGeometry))
{
// now, create entries for all fracture scvs on this same facet ...
if (scv.isOnFracture() && scv.facetIndexInElement() == thisFacetIdx)
jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
// ... and those non-fracture scvs that are not on this facet
else if (!scv.isOnFracture()
&& std::find( facetLocalDofs.begin(),
facetLocalDofs.end(),
scv.localDofIndex() ) == facetLocalDofs.end())
{
jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
}
}
}
}
//! return the Jacobian of the shape functions at the integration point
const std::vector<ShapeJacobian>& shapeJacobian() const { return shapeJacobian_; }
//! return the shape values for all scvs at the integration point
const std::vector<ShapeValue>& shapeValues() const { return shapeValues_; }
//! return the shape value gradients for all scvs at the integration point
const JacobianInverseTransposed& jacInvT() const { return jacInvT_; }
//! return the shape value gradients corresponding to an scv
const GlobalPosition& gradN(unsigned int scvIdxInElement) const { return gradN_[scvIdxInElement]; }
private:
std::vector<GlobalPosition> gradN_;
std::vector<ShapeJacobian> shapeJacobian_;
std::vector<ShapeValue> shapeValues_;
JacobianInverseTransposed jacInvT_;
};
} // end namespace Dumux
#endif
This diff is collapsed.
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 2 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 BoxDiscretization
* \ingroup BoxDFM
* \brief Helper class constructing the dual grid finite volume geometries
* for the box discrete fracture model.
*/
#ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_GEOMETRY_HELPER_HH
#define DUMUX_POROUSMEDIUMFLOW_BOXDFM_GEOMETRY_HELPER_HH
#include <dumux/discretization/box/boxgeometryhelper.hh>
namespace Dumux
{
//! Create sub control volumes and sub control volume face geometries
template<class GridView, int dim, class ScvType, class ScvfType>
class BoxDfmGeometryHelper;
//! A class to create sub control volume and sub control volume face geometries per element
template <class GridView, class ScvType, class ScvfType>
class BoxDfmGeometryHelper<GridView, 2, ScvType, ScvfType> : public BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
{
using ParentType = BoxGeometryHelper<GridView, 2, ScvType, ScvfType>;
using Intersection = typename GridView::Intersection;
using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
static constexpr auto dim = GridView::dimension;
using Scalar = typename GridView::ctype;
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
public:
//! Pull up constructor of base class
using ParentType::ParentType;
//! Get the corners of the (d-1)-dimensional fracture scvf
//! The second argument is for compatibility reasons with the 3d case!
ScvfCornerStorage getFractureScvfCorners(const Intersection& is,
const typename Intersection::Geometry& isGeom,
unsigned int idxOnIntersection = 0) const
{ return ScvfCornerStorage({isGeom.center()}); }
//! get fracture scvf normal vector (simply the unit vector of the edge)
//! The third argument is for compatibility reasons with the 3d case!
typename ScvfType::Traits::GlobalPosition
fractureNormal(const ScvfCornerStorage& p,
const Intersection& is,
unsigned int edgeIndexInIntersection = 0) const
{
assert(isGeom.corners() == 2);
const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
const auto vIdxLocal0 = referenceElement.subEntity(is.indexInInside(), 1, 0, dim);
const auto vIdxLocal1 = referenceElement.subEntity(is.indexInInside(), 1, 1, dim);
auto n = this->elementGeometry_.corner(vIdxLocal1) - this->elementGeometry_.corner(vIdxLocal0);
n /= n.two_norm();
return n;
}
};
//! A class to create sub control volume and sub control volume face geometries per element
template <class GridView, class ScvType, class ScvfType>
class BoxDfmGeometryHelper<GridView, 3, ScvType, ScvfType> : public BoxGeometryHelper<GridView, 3, ScvType, ScvfType>
{
using ParentType = BoxGeometryHelper<GridView, 3, ScvType, ScvfType>;
using Intersection = typename GridView::Intersection;
using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
static constexpr auto dim = GridView::dimension;
static constexpr auto dimWorld = GridView::dimensionworld;
using Scalar = typename GridView::ctype;
using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
public:
//! Pull up constructor of base class
using ParentType::ParentType;
//! Create the sub control volume face geometries on an intersection marked as fracture
ScvfCornerStorage getFractureScvfCorners(const Intersection& is,
const typename Intersection::Geometry& isGeom,
unsigned int edgeIndexInIntersection) const
{
const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
const auto faceRefElem = FaceReferenceElements::general(isGeom.type());
// create point vector for this geometry
typename ScvfType::Traits::GlobalPosition pi[9];
// the facet center
pi[0] = isGeom.center();
// facet edge midpoints
const auto idxInInside = is.indexInInside();
for (int i = 0; i < faceRefElem.size(1); ++i)
{
const auto edgeIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim-1);
pi[i+1] = this->p[edgeIdxLocal+this->corners_+1];
}
// proceed according to number of corners
const auto corners = isGeom.corners();
switch (corners)
{
case 3: // triangle
{
//! Only build the maps the first time we encounter a triangle
static const std::uint8_t fo = 1; //!< face offset in point vector p
static const std::uint8_t map[3][2] =
{
{fo+0, 0},
{0, fo+1},
{fo+2, 0}
};
return ScvfCornerStorage{ {pi[map[edgeIndexInIntersection][0]],
pi[map[edgeIndexInIntersection][1]]} };
}
case 4: // quadrilateral
{
//! Only build the maps the first time we encounter a quadrilateral
static const std::uint8_t fo = 1; //!< face offset in point vector p
static const std::uint8_t map[4][2] =
{
{0, fo+0},
{fo+1, 0},
{fo+2, 0},
{0, fo+3}
};
return ScvfCornerStorage{ {pi[map[edgeIndexInIntersection][0]],
pi[map[edgeIndexInIntersection][1]]} };
}
default:
DUNE_THROW(Dune::NotImplemented, "Box fracture scvf geometries for dim=" << dim
<< " dimWorld=" << dimWorld
<< " corners=" << corners);
}
}
//! get fracture scvf normal vector
typename ScvfType::Traits::GlobalPosition
fractureNormal(const ScvfCornerStorage& p,
const Intersection& is,
unsigned int edgeIndexInIntersection) const
{
const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
// first get the intersection corners (maximum "4" is for quadrilateral face)
typename ScvfType::Traits::GlobalPosition c[4];
const auto corners = referenceElement.size(is.indexInInside(), 1, dim);
for (int i = 0; i < corners; ++i)
{
const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, i, dim);
c[i] = this->elementGeometry_.corner(vIdxLocal);
}
// compute edge vector depending on number of corners
const auto gridEdge = [&] ()
{
// triangles
if (corners == 3)
{
if (edgeIndexInIntersection == 0) return c[1]-c[0];
else if (edgeIndexInIntersection == 1) return c[2]-c[0];
else if (edgeIndexInIntersection == 2) return c[2]-c[1];
else DUNE_THROW(Dune::InvalidStateException, "Invalid edge index");
}
else if (corners == 4)
{
if (edgeIndexInIntersection == 0) return c[2]-c[0];
else if (edgeIndexInIntersection == 1) return c[3]-c[1];
else if (edgeIndexInIntersection == 2) return c[1]-c[0];
else if (edgeIndexInIntersection == 3) return c[3]-c[2];
else DUNE_THROW(Dune::InvalidStateException, "Invalid edge index");
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid face geometry");
} ();
// compute lower edge of the scvf
assert(p.size() == 2);
const auto scvfEdge = p[1]-p[0];
// compute scvf normal via 2 cross products
const auto faceN = crossProduct(gridEdge, scvfEdge);
auto n = crossProduct(scvfEdge, faceN);
n /= n.two_norm();
return n;
}
};
} // end namespace Dumux
#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 2 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/>. *
*****************************************************************************/
/*!
* \ingroup Properties
* \file
* \brief Defines a type tag and some properties for porous medium
* flow models using the box scheme extended to discrete fractures.
*/
#ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_MODEL_HH
#define DUMUX_POROUSMEDIUMFLOW_BOXDFM_MODEL_HH
#include <dumux/discretization/box/properties.hh>
#include "fvgridgeometry.hh"
#include "fluxvariablescache.hh"
namespace Dumux {
namespace Properties {
//! Type tag for the box scheme.
NEW_TYPE_TAG(BoxDfmModel, INHERITS_FROM(BoxModel));
//! Set the default for the global finite volume geometry
SET_PROP(BoxDfmModel, FVGridGeometry)
{
private:
static constexpr bool enableCache = GET_PROP_VALUE(TypeTag, EnableFVGridGeometryCache);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = BoxDfmFVGridGeometry<Scalar, GridView, enableCache>;
};
//! The flux variables cache class specific to box-dfm porous medium flow models
SET_TYPE_PROP(BoxDfmModel, FluxVariablesCache, BoxDfmFluxVariablesCache<TypeTag>);
} // namespace Properties
} // namespace Dumux
#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 2 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 BoxDiscretization
* \brief the sub control volume for the box discrete fracture scheme
*/
#ifndef DUMUX_POROUSMEDIUMFLOW_BOXDFM_SUBCONTROLVOLUME_HH
#define DUMUX_POROUSMEDIUMFLOW_BOXDFM_SUBCONTROLVOLUME_HH
#include <dune/common/version.hh>
#include <dune/common/reservedvector.hh>
#include <dumux/discretization/subcontrolvolumebase.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
#include <dumux/common/math.hh>
namespace Dumux {
/*!
* \ingroup BoxDiscretization
* \brief Default traits class to be used for the sub-control volumes
* for the box discrete fracture scheme
* \tparam GV the type of the grid view
*
* \note We define new traits for the box-dfm sub-control volume face
* as we use a different type of container for storing the scvf corners!
*/
template<class GridView>
struct BoxDfmDefaultScvGeometryTraits
{
using Grid = typename GridView::Grid;
static const int dim = Grid::dimension;
static const int dimWorld = Grid::dimensionworld;
template <class ct>
struct ScvMLGTraits : public Dune::MultiLinearGeometryTraits<ct>
{
// we use static vectors to store the corners as we know
// the number of corners in advance (2^(dim) corners (1<<(dim))
// However, on fracture scvs the number might be smaller (use ReservedVector)
template< int mydim, int cdim >
struct CornerStorage
{
using Type = Dune::ReservedVector< Dune::FieldVector< ct, cdim >, (1<<(dim)) >;
};
// we know all scvfs will have the same geometry type
template< int mydim >
struct hasSingleGeometryType
{
static const bool v = true;
static const unsigned int topologyId = Dune::Impl::CubeTopology< mydim >::type::id;
};
};
using GridIndexType = typename Grid::LeafGridView::IndexSet::IndexType;
using LocalIndexType = unsigned int;
using Scalar = typename Grid::ctype;
using Geometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld, ScvMLGTraits<Scalar>>;
using CornerStorage = typename ScvMLGTraits<Scalar>::template CornerStorage<dim, dimWorld>::Type;
using GlobalPosition = typename CornerStorage::value_type;
};
/*!
* \ingroup BoxDiscretization
* \brief the sub control volume for the box discrete fracture scheme
* \tparam GV the type of the grid view
* \tparam T the scvf geometry traits
*/
template<class GV,
class T = BoxDfmDefaultScvGeometryTraits<GV> >
class BoxDfmSubControlVolume
: public SubControlVolumeBase<BoxDfmSubControlVolume<GV, T>, T>
{
using ThisType = BoxDfmSubControlVolume<GV, T>;
using ParentType = SubControlVolumeBase<ThisType, T>;
using Geometry = typename T::Geometry;
using GridIndexType = typename T::GridIndexType;
using LocalIndexType = typename T::LocalIndexType;
using Scalar = typename T::Scalar;
using GlobalPosition = typename T::GlobalPosition;
using CornerStorage = typename T::CornerStorage;
enum { dim = Geometry::mydimension };
static_assert(dim == 2 || dim == 3, "Box-Dfm sub-control volume only implemented in 2d or 3d");
public:
//! state the traits public and thus export all types
using Traits = T;