Merge branch 'feature/boxdfm' into 'master'

Feature/boxdfm

See merge request !946

dumux/discretization/box/boxgeometryhelper.hh

@@ -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

dumux/io/vtkoutputmodule.hh

@@ -235,6 +235,24 @@ public:
         }
     }
 
+protected:
+    // some return functions for differing implementations to use
+    const Problem& problem() const { return problem_; }
+    const FVGridGeometry& fvGridGeometry() const { return gridGeom_; }
+    const GridVariables& gridVariables() const { return gridVariables_; }
+    const SolutionVector& sol() const { return sol_; }
+
+    const bool verbose() const { return verbose_; }
+    const std::string& name() const { return name_; }
+    const Dune::VTK::DataMode dataMode() const { return dm_; }
+
+    Dune::VTKWriter<GridView>& writer() { return *writer_; }
+    Dune::VTKSequenceWriter<GridView>& sequenceWriter() { return sequenceWriter_; }
+
+    const std::vector<VolVarScalarDataInfo>& volVarScalarDataInfo() const { return volVarScalarDataInfo_; }
+    const std::vector<VolVarVectorDataInfo>& volVarVectorDataInfo() const { return volVarVectorDataInfo_; }
+    const std::vector<Field>& fields() const { return fields_; }
+
 private:
 
     //! Assembles the fields and adds them to the writer (conforming output)

dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh

@@ -26,6 +26,8 @@
 #ifndef DUMUX_BROOKS_COREY_PARAMS_HH
 #define DUMUX_BROOKS_COREY_PARAMS_HH
 
+#include <dune/common/float_cmp.hh>
+
 #include <dumux/common/valgrind.hh>
 
 namespace Dumux
@@ -55,6 +57,16 @@ public:
     {
     }
 
+    /*!
+     * \brief Equality comparison with another set of params
+     */
+    template<class OtherParams>
+    bool operator== (const OtherParams& otherParams) const
+    {
+        return Dune::FloatCmp::eq(pe_, otherParams.pe(), /*eps*/1e-6*pe_)
+               && Dune::FloatCmp::eq(lambda_, otherParams.lambda(), /*eps*/1e-6*lambda_);
+    }
+
     /*!
      * \brief Returns the entry pressure in \f$\mathrm{[Pa]}\f$
      */

dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh

@@ -25,6 +25,8 @@
 #ifndef DUMUX_REGULARIZED_BROOKS_COREY_PARAMS_HH
 #define DUMUX_REGULARIZED_BROOKS_COREY_PARAMS_HH
 
+#include <dune/common/float_cmp.hh>
+
 #include "brookscoreyparams.hh"
 
 namespace Dumux
@@ -54,6 +56,16 @@ public:
         setThresholdSw(0.01);
     }
 
+    /*!
+     * \brief Equality comparison with another set of params
+     */
+    template<class OtherParams>
+    bool operator== (const OtherParams& otherParams) const
+    {
+        return Dune::FloatCmp::eq(thresholdSw_, otherParams.thresholdSw(), /*eps*/1e-6*thresholdSw_)
+               && BrooksCoreyParams::operator==(otherParams);
+    }
+
     /*!
      * \brief Set the threshold saturation below which the capillary pressure
      *        is regularized.

dumux/porousmediumflow/CMakeLists.txt

@@ -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")

dumux/porousmediumflow/boxdfm/CMakeLists.txt

+
+#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)

dumux/porousmediumflow/boxdfm/fluxvariablescache.hh

+// -*- 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

dumux/porousmediumflow/boxdfm/geometryhelper.hh

+// -*- 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

dumux/porousmediumflow/boxdfm/model.hh

+// -*- 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

dumux/porousmediumflow/boxdfm/subcontrolvolume.hh

+// -*- 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