diff --git a/dumux/implicit/cornerpoint/cpdarcyfluxvariables.hh b/dumux/implicit/cornerpoint/cpdarcyfluxvariables.hh
index 7a56d40abef95467e7e2a06588cd8599457dea31..06dd4e0243c6019f531e4883e3f32e0ddd30e909 100644
--- a/dumux/implicit/cornerpoint/cpdarcyfluxvariables.hh
+++ b/dumux/implicit/cornerpoint/cpdarcyfluxvariables.hh
@@ -1,278 +1,8 @@
-// -*- 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 This file contains the data which is required to calculate
- * volume fluxes of fluid phases over a face of a finite volume by means
- * of the Darcy approximation.
- *
- */
-#ifndef DUMUX_CP_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_CP_DARCY_FLUX_VARIABLES_HH
+#ifndef DUMUX_CP_DARCY_FLUX_VARIABLES_HH_OLD
+#define DUMUX_CP_DARCY_FLUX_VARIABLES_HH_OLD
-#include <dune/common/float_cmp.hh>
+#warning this header is deprecated, use dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh instead
-#include <dumux/common/math.hh>
-#include <dumux/common/parameters.hh>
+#include <dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh>
-#include <dumux/implicit/properties.hh>
-
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}
-
-/*!
- * \ingroup ImplicitFluxVariables
- * \brief Evaluates the normal component of the Darcy velocity
- * on a (sub)control volume face.
- */
-template <class TypeTag>
-class CpDarcyFluxVariables
-{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
- enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldVector<Scalar, dim> DimVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*!
- * \brief The constructor
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- CpDarcyFluxVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- : fvGeometry_(fvGeometry), faceIdx_(fIdx), onBoundary_(onBoundary)
- {
- mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
- calculateVolumeFlux_(problem, element, elemVolVars);
- }
-
-public:
- /*!
- * \brief Return the volumetric flux over a face of a given phase.
- *
- * This is the calculated velocity multiplied by the unit normal
- * and the area of the face. face().normal has already the
- * magnitude of the area.
- *
- * \param phaseIdx index of the phase
- */
- Scalar volumeFlux(const unsigned int phaseIdx) const
- { return volumeFlux_[phaseIdx]; }
-
- /*!
- * \brief Return the velocity of a given phase.
- *
- * This is the full velocity vector on the
- * face (without being multiplied with normal).
- *
- * \param phaseIdx index of the phase
- */
- GlobalPosition velocity(const unsigned int phaseIdx) const
- { return velocity_[phaseIdx] ; }
-
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- const unsigned int downstreamIdx(const unsigned phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
- *
- * \param phaseIdx index of the phase
- */
- const unsigned int upstreamIdx(const unsigned phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
-
- /*!
- * \brief Return the SCV (sub-control-volume) face. This may be either
- * a face within the element or a face on the element boundary,
- * depending on the value of onBoundary_.
- */
- const SCVFace &face() const
- {
- if (onBoundary_)
- return fvGeometry_.boundaryFace[faceIdx_];
- else
- return fvGeometry_.subContVolFace[faceIdx_];
- }
-
-protected:
- /*!
- * \brief Actual calculation of the volume flux.
- *
- * \param problem The problem
- * \param element The finite element
- * \param elemVolVars The volume variables of the current element
- */
- void calculateVolumeFlux_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate the transmissibilities
- const SpatialParams &spatialParams = problem.spatialParams();
-
- const Element& elementI = fvGeometry_.neighbors[face().i];
- FVElementGeometry fvGeometryI;
- fvGeometryI.subContVol[0].global = elementI.geometry().center();
- auto ki = spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0);
- Dune::FieldVector<Scalar, dimWorld> kin;
- ki.mv(face().normal, kin);
- kin /= face().area;
- auto di = face().ipGlobal;
- di -= elementI.geometry().center();
- auto ti = std::abs(di*kin*face().area/(2*di.two_norm2()));
-
- auto tij = ti;
- if (!onBoundary_)
- {
- const Element& elementJ = fvGeometry_.neighbors[face().j];
- FVElementGeometry fvGeometryJ;
- fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
- auto kj = spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0);
- Dune::FieldVector<Scalar, dimWorld> kjn;
- kj.mv(face().normal, kjn);
- kjn /= face().area;
- auto dj = face().ipGlobal;
- dj -= elementJ.geometry().center();
- auto tj = std::abs(dj*kjn*face().area/(2*dj.two_norm2()));
- tij = harmonicMean(ti, tj);
- }
-
- // loop over all phases
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- const auto& volVarsI = elemVolVars[face().i];
- auto potentialI = volVarsI.pressure(phaseIdx);
-
- const auto& volVarsJ = elemVolVars[face().j];
- auto potentialJ = volVarsJ.pressure(phaseIdx);
-
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
- {
- // calculate the phase density at the integration point. we
- // only do this if the phase is present in both cells
- Scalar SI = volVarsI.fluidState().saturation(phaseIdx);
- Scalar SJ = volVarsJ.fluidState().saturation(phaseIdx);
- Scalar rhoI = volVarsI.fluidState().density(phaseIdx);
- Scalar rhoJ = volVarsJ.fluidState().density(phaseIdx);
- Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
- Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
- if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
- // doesn't matter because no wetting phase is present in
- // both cells!
- fI = fJ = 0.5;
- const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
-
- auto globalPosI = elementI.geometry().center();
- potentialI -= density*(problem.gravityAtPos(globalPosI)*globalPosI);
-
- if (onBoundary_)
- {
- potentialJ -= density*(problem.gravityAtPos(face().ipGlobal)*face().ipGlobal);
- }
- else
- {
- const Element& elementJ = fvGeometry_.neighbors[face().j];
- auto globalPosJ = elementJ.geometry().center();
- potentialJ -= density*(problem.gravityAtPos(globalPosJ)*globalPosJ);
- }
- }
-
- auto potentialDiff = potentialI - potentialJ;
-
- // determine the upwind direction
- if (potentialDiff > 0)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
- }
- else
- {
- upstreamIdx_[phaseIdx] = face().j;
- downstreamIdx_[phaseIdx] = face().i;
- }
-
- // obtain the upwind volume variables
- const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
- const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
-
- // set the volume flux
- volumeFlux_[phaseIdx] = tij*potentialDiff;
- volumeFlux_[phaseIdx] *= mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
- + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx);
-
- velocity_[phaseIdx] = face().normal;
- velocity_[phaseIdx] *= volumeFlux_[phaseIdx]/face().area;
- } // over loop all phases
- }
-
- const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
- const unsigned int faceIdx_; //!< The index of the sub control volume face
- const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
- unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
- Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
- GlobalPosition velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
- Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
-};
-
-} // end namespace
-
-#endif // DUMUX_CP_DARCY_FLUX_VARIABLES_HH
+#endif
diff --git a/dumux/implicit/cornerpoint/cpelementvolumevariables.hh b/dumux/implicit/cornerpoint/cpelementvolumevariables.hh
index 7a8da16acb31461fb644d89caf2d220d61d69997..a2e13dfb105a7808653304e9cf0e81719dbc747f 100644
--- a/dumux/implicit/cornerpoint/cpelementvolumevariables.hh
+++ b/dumux/implicit/cornerpoint/cpelementvolumevariables.hh
@@ -1,142 +1,8 @@
-// -*- 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 Volume variables gathered on an element
- */
-#ifndef DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
-#define DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
+#ifndef DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH_OLD
+#define DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH_OLD
-#include <dumux/implicit/cellcentered/properties.hh>
+#warning this header is deprecated, use dumux/implicit/cornerpoint/elementvolumevariables.hh instead
-namespace Dumux
-{
-
-/*!
- * \ingroup CCModel
- * \brief This class stores an array of VolumeVariables objects, one
- * volume variables object for each of the element's vertices
- */
-template<class TypeTag>
-class CpElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) >
-{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
-public:
- /*!
- * \brief The constructor.
- */
- CpElementVolumeVariables()
- { }
-
- /*!
- * \brief Construct the volume variables for all of vertices of an element.
- *
- * \param problem The problem which needs to be simulated.
- * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvGeometry The finite volume geometry of the element
- * \param oldSol Tells whether the model's previous or current solution should be used.
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- bool oldSol)
- {
- const SolutionVector &globalSol =
- oldSol?
- problem.model().prevSol():
- problem.model().curSol();
-
- int numNeighbors = fvGeometry.numNeighbors;
- this->resize(numNeighbors);
-
- for (int i = 0; i < numNeighbors; i++)
- {
- const Element& neighbor = fvGeometry.neighbors[i];
- const PrimaryVariables &solI
- = globalSol[problem.elementMapper().index(neighbor)];
-
- FVElementGeometry neighborFVGeom;
- neighborFVGeom.updateInner(neighbor);
-
- (*this)[i].update(solI,
- problem,
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- oldSol);
- }
-
- // only treat boundary if current solution is evaluated
- if (!oldSol)
- {
- // check if element intersects with the boundary
- ElementBoundaryTypes elemBCTypes;
- elemBCTypes.update(problem, element);
- if (elemBCTypes.hasDirichlet()
- || elemBCTypes.hasNeumann()
- || elemBCTypes.hasOutflow())
- {
- const int numFaces = 6;
- this->resize(numNeighbors + numFaces);
-
- // add volume variables for the boundary faces
- for (const auto& intersection : Dune::intersections(problem.gridView(), element)) {
- if (!intersection.boundary())
- continue;
-
- BoundaryTypes bcTypes;
- problem.boundaryTypes(bcTypes, intersection);
-
- int fIdx = intersection.indexInInside();
- int indexInVariables = numNeighbors + fIdx;
-
- if (bcTypes.hasDirichlet())
- {
- PrimaryVariables dirichletValues;
- problem.dirichlet(dirichletValues, intersection);
-
- (*this)[indexInVariables].update(dirichletValues,
- problem,
- element,
- fvGeometry,
- /*scvIdx=*/0,
- oldSol);
- }
- else
- {
- (*this)[indexInVariables] = (*this)[0];
- }
- }
- }
- }
- }
-};
-
-} // namespace Dumux
+#include <dumux/implicit/cornerpoint/elementvolumevariables.hh>
#endif
diff --git a/dumux/implicit/cornerpoint/cpfvelementgeometry.hh b/dumux/implicit/cornerpoint/cpfvelementgeometry.hh
index d60d0604a883d5b82b79cd371547e9f482732a3a..c35ef0758db9d033b5c1e165be1e5100ca22e821 100644
--- a/dumux/implicit/cornerpoint/cpfvelementgeometry.hh
+++ b/dumux/implicit/cornerpoint/cpfvelementgeometry.hh
@@ -1,233 +1,8 @@
-// -*- 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 Represents the finite volume geometry of a single element in
- * the cell-centered fv scheme.
- */
-#ifndef DUMUX_CP_FV_ELEMENTGEOMETRY_HH
-#define DUMUX_CP_FV_ELEMENTGEOMETRY_HH
+#ifndef DUMUX_CP_FV_ELEMENTGEOMETRY_HH_OLD
+#define DUMUX_CP_FV_ELEMENTGEOMETRY_HH_OLD
-#include <dune/geometry/referenceelements.hh>
-#include <dune/grid/common/intersectioniterator.hh>
+#warning this header is deprecated, use dumux/implicit/cornerpoint/fvelementgeometry.hh instead
-#include <dumux/common/propertysystem.hh>
-
-namespace Dumux
-{
-namespace Properties
-{
-NEW_PROP_TAG(GridView);
-NEW_PROP_TAG(Scalar);
-}
-
-/*!
- * \ingroup CCModel
- * \brief Represents the finite volume geometry of a single element in
- * the cell-centered fv scheme.
- */
-template<class TypeTag>
-class CpFVElementGeometry
-{
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- enum{dim = GridView::dimension};
- enum{dimWorld = GridView::dimensionworld};
-
- enum{maxNFAP = 2}; //! maximum number of flux approximation points (two-point flux)
- enum{maxNE = 50}; //! maximum number of neighbors
- enum{maxBF = 2*dim}; //! maximum number of boundary faces
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GridView::ctype CoordScalar;
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename Element::Geometry Geometry;
- typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
- typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
-
-public:
- struct SubControlVolume //! FV intersected with element
- {
- LocalPosition local; //!< local position
- GlobalPosition global; //!< global position
- Scalar volume; //!< volume of scv
- bool inner;
- };
-
- struct SubControlVolumeFace //! interior face of a sub control volume
- {
- int i,j; //!< scvf seperates corner i and j of elem
- LocalPosition ipLocal; //!< integration point in local coords
- GlobalPosition ipGlobal; //!< integration point in global coords
- GlobalPosition normal; //!< normal on face pointing to CV j or outward of the domain with length equal to |scvf|
- Scalar area; //!< area of face
- Dune::FieldVector<GlobalPosition, maxNFAP> grad; //!< derivatives of shape functions at ip
- Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip
- Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points
- unsigned numFap; //!< number of flux approximation points
- unsigned fIdx; //!< the index (w.r.t. the element) of the face (codim 1 entity) that the scvf is part of
- };
-
- typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef
-
- LocalPosition elementLocal; //!< local coordinate of element center
- GlobalPosition elementGlobal; //!< global coordinate of element center
- Scalar elementVolume; //!< element volume
- SubControlVolume subContVol[1]; //!< data of the sub control volumes
- SubControlVolumeFace subContVolFace[maxNE]; //!< data of the sub control volume faces
- BoundaryFace boundaryFace[maxBF]; //!< data of the boundary faces
- int numScv; //!< number of subcontrol volumes
- int numScvf; //!< number of inner-domain subcontrolvolume faces
- int numNeighbors; //!< number of neighboring elements including the element itself
- std::vector<Element> neighbors; //!< stores the neighboring elements
-
- void updateInner(const Element& element)
- {
- const Geometry geometry = element.geometry();
-
- elementVolume = geometry.volume();
- elementGlobal = geometry.center();
- //elementLocal = geometry.local(elementGlobal);
-
- numScv = 1;
- numScvf = 0;
-
- subContVol[0].local = elementLocal;
- subContVol[0].global = elementGlobal;
- subContVol[0].inner = true;
- subContVol[0].volume = elementVolume;
-
- // initialize neighbors list with self:
- numNeighbors = 1;
- neighbors.clear();
- neighbors.reserve(maxNE);
- neighbors.push_back(element);
- }
-
- void update(const GridView& gridView, const Element& element)
- {
- updateInner(element);
-
- const Geometry geometry = element.geometry();
-
- bool onBoundary = false;
-
- // fill neighbor information and control volume face data:
- for (const auto& intersection : Dune::intersections(gridView, element))
- {
- const auto isGeometry = intersection.geometry();
-
- // neighbor information and inner cvf data:
- if (intersection.neighbor())
- {
- numNeighbors++;
- neighbors.push_back(intersection.outside());
-
- int scvfIdx = numNeighbors - 2;
- SubControlVolumeFace& scvFace = subContVolFace[scvfIdx];
-
- scvFace.i = 0;
- scvFace.j = scvfIdx + 1;
-
- scvFace.ipGlobal = isGeometry.center();
- //scvFace.ipLocal = geometry.local(scvFace.ipGlobal);
- scvFace.normal = intersection.centerUnitOuterNormal();
- auto di = scvFace.ipGlobal;
- di -= elementGlobal;
- if (scvFace.normal*di < 0)
- scvFace.normal *= -1.0;
- Scalar volume = isGeometry.volume();
- scvFace.normal *= volume;
- scvFace.area = volume;
-
- GlobalPosition distVec = elementGlobal
- - neighbors[scvfIdx+1].geometry().center();
- distVec /= distVec.two_norm2();
-
- // gradients using a two-point flux approximation
- scvFace.numFap = 2;
- for (unsigned int fapIdx = 0; fapIdx < scvFace.numFap; fapIdx++)
- {
- scvFace.grad[fapIdx] = distVec;
- scvFace.shapeValue[fapIdx] = 0.5;
- }
- scvFace.grad[1] *= -1.0;
-
- scvFace.fapIndices[0] = scvFace.i;
- scvFace.fapIndices[1] = scvFace.j;
-
- scvFace.fIdx = intersection.indexInInside();
- }
-
- // boundary cvf data
- if (intersection.boundary())
- {
- onBoundary = true;
- int bfIdx = intersection.indexInInside();
- SubControlVolumeFace& bFace = boundaryFace[bfIdx];
-
- bFace.ipGlobal = isGeometry.center();
- //bFace.ipLocal = geometry.local(bFace.ipGlobal);
- bFace.normal = intersection.centerUnitOuterNormal();
- auto di = bFace.ipGlobal;
- di -= elementGlobal;
- if (bFace.normal*di < 0)
- bFace.normal *= -1.0;
- Scalar volume = isGeometry.volume();
- bFace.normal *= volume;
- bFace.area = volume;
- bFace.i = 0;
- bFace.j = 0;
-
- GlobalPosition distVec = elementGlobal - bFace.ipGlobal;
- distVec /= distVec.two_norm2();
-
- // gradients using a two-point flux approximation
- bFace.numFap = 2;
- for (unsigned int fapIdx = 0; fapIdx < bFace.numFap; fapIdx++)
- {
- bFace.grad[fapIdx] = distVec;
- bFace.shapeValue[fapIdx] = 0.5;
- }
- bFace.grad[1] *= -1.0;
-
- bFace.fapIndices[0] = bFace.i;
- bFace.fapIndices[1] = bFace.j;
- }
- }
-
- // set the number of inner-domain subcontrolvolume faces
- numScvf = numNeighbors - 1;
-
- // treat elements on the boundary
- if (onBoundary)
- {
- for (int bfIdx = 0; bfIdx < maxBF; bfIdx++)
- {
- SubControlVolumeFace& bFace = boundaryFace[bfIdx];
- bFace.j = numNeighbors + bfIdx;
- bFace.fapIndices[1] = bFace.j;
- neighbors.push_back(element);
- }
- }
- }
-};
-
-}
+#include <dumux/implicit/cornerpoint/fvelementgeometry.hh>
#endif
-
diff --git a/dumux/implicit/cornerpoint/elementvolumevariables.hh b/dumux/implicit/cornerpoint/elementvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7a8da16acb31461fb644d89caf2d220d61d69997
--- /dev/null
+++ b/dumux/implicit/cornerpoint/elementvolumevariables.hh
@@ -0,0 +1,142 @@
+// -*- 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 Volume variables gathered on an element
+ */
+#ifndef DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
+#define DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
+
+#include <dumux/implicit/cellcentered/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup CCModel
+ * \brief This class stores an array of VolumeVariables objects, one
+ * volume variables object for each of the element's vertices
+ */
+template<class TypeTag>
+class CpElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) >
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+public:
+ /*!
+ * \brief The constructor.
+ */
+ CpElementVolumeVariables()
+ { }
+
+ /*!
+ * \brief Construct the volume variables for all of vertices of an element.
+ *
+ * \param problem The problem which needs to be simulated.
+ * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
+ * \param fvGeometry The finite volume geometry of the element
+ * \param oldSol Tells whether the model's previous or current solution should be used.
+ */
+ void update(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ bool oldSol)
+ {
+ const SolutionVector &globalSol =
+ oldSol?
+ problem.model().prevSol():
+ problem.model().curSol();
+
+ int numNeighbors = fvGeometry.numNeighbors;
+ this->resize(numNeighbors);
+
+ for (int i = 0; i < numNeighbors; i++)
+ {
+ const Element& neighbor = fvGeometry.neighbors[i];
+ const PrimaryVariables &solI
+ = globalSol[problem.elementMapper().index(neighbor)];
+
+ FVElementGeometry neighborFVGeom;
+ neighborFVGeom.updateInner(neighbor);
+
+ (*this)[i].update(solI,
+ problem,
+ neighbor,
+ neighborFVGeom,
+ /*scvIdx=*/0,
+ oldSol);
+ }
+
+ // only treat boundary if current solution is evaluated
+ if (!oldSol)
+ {
+ // check if element intersects with the boundary
+ ElementBoundaryTypes elemBCTypes;
+ elemBCTypes.update(problem, element);
+ if (elemBCTypes.hasDirichlet()
+ || elemBCTypes.hasNeumann()
+ || elemBCTypes.hasOutflow())
+ {
+ const int numFaces = 6;
+ this->resize(numNeighbors + numFaces);
+
+ // add volume variables for the boundary faces
+ for (const auto& intersection : Dune::intersections(problem.gridView(), element)) {
+ if (!intersection.boundary())
+ continue;
+
+ BoundaryTypes bcTypes;
+ problem.boundaryTypes(bcTypes, intersection);
+
+ int fIdx = intersection.indexInInside();
+ int indexInVariables = numNeighbors + fIdx;
+
+ if (bcTypes.hasDirichlet())
+ {
+ PrimaryVariables dirichletValues;
+ problem.dirichlet(dirichletValues, intersection);
+
+ (*this)[indexInVariables].update(dirichletValues,
+ problem,
+ element,
+ fvGeometry,
+ /*scvIdx=*/0,
+ oldSol);
+ }
+ else
+ {
+ (*this)[indexInVariables] = (*this)[0];
+ }
+ }
+ }
+ }
+ }
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/implicit/cornerpoint/fvelementgeometry.hh b/dumux/implicit/cornerpoint/fvelementgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d25af5ec8182f763584cca4416d28013939dfdf9
--- /dev/null
+++ b/dumux/implicit/cornerpoint/fvelementgeometry.hh
@@ -0,0 +1,232 @@
+// -*- 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 Represents the finite volume geometry of a single element in
+ * the cell-centered fv scheme.
+ */
+#ifndef DUMUX_CP_FV_ELEMENTGEOMETRY_HH
+#define DUMUX_CP_FV_ELEMENTGEOMETRY_HH
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/grid/common/intersectioniterator.hh>
+
+#include <dumux/common/propertysystem.hh>
+
+namespace Dumux
+{
+namespace Properties
+{
+NEW_PROP_TAG(GridView);
+NEW_PROP_TAG(Scalar);
+}
+
+/*!
+ * \ingroup CCModel
+ * \brief Represents the finite volume geometry of a single element in
+ * the cell-centered fv scheme.
+ */
+template<class TypeTag>
+class CpFVElementGeometry
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ enum{dim = GridView::dimension};
+ enum{dimWorld = GridView::dimensionworld};
+
+ enum{maxNFAP = 2}; //! maximum number of flux approximation points (two-point flux)
+ enum{maxNE = 50}; //! maximum number of neighbors
+ enum{maxBF = 2*dim}; //! maximum number of boundary faces
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename Element::Geometry Geometry;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+
+public:
+ struct SubControlVolume //! FV intersected with element
+ {
+ LocalPosition local; //!< local position
+ GlobalPosition global; //!< global position
+ Scalar volume; //!< volume of scv
+ bool inner;
+ };
+
+ struct SubControlVolumeFace //! interior face of a sub control volume
+ {
+ int i,j; //!< scvf seperates corner i and j of elem
+ LocalPosition ipLocal; //!< integration point in local coords
+ GlobalPosition ipGlobal; //!< integration point in global coords
+ GlobalPosition normal; //!< normal on face pointing to CV j or outward of the domain with length equal to |scvf|
+ Scalar area; //!< area of face
+ Dune::FieldVector<GlobalPosition, maxNFAP> grad; //!< derivatives of shape functions at ip
+ Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip
+ Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points
+ unsigned numFap; //!< number of flux approximation points
+ unsigned fIdx; //!< the index (w.r.t. the element) of the face (codim 1 entity) that the scvf is part of
+ };
+
+ typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef
+
+ LocalPosition elementLocal; //!< local coordinate of element center
+ GlobalPosition elementGlobal; //!< global coordinate of element center
+ Scalar elementVolume; //!< element volume
+ SubControlVolume subContVol[1]; //!< data of the sub control volumes
+ SubControlVolumeFace subContVolFace[maxNE]; //!< data of the sub control volume faces
+ BoundaryFace boundaryFace[maxBF]; //!< data of the boundary faces
+ int numScv; //!< number of subcontrol volumes
+ int numScvf; //!< number of inner-domain subcontrolvolume faces
+ int numNeighbors; //!< number of neighboring elements including the element itself
+ std::vector<Element> neighbors; //!< stores the neighboring elements
+
+ void updateInner(const Element& element)
+ {
+ const Geometry geometry = element.geometry();
+
+ elementVolume = geometry.volume();
+ elementGlobal = geometry.center();
+ //elementLocal = geometry.local(elementGlobal);
+
+ numScv = 1;
+ numScvf = 0;
+
+ subContVol[0].local = elementLocal;
+ subContVol[0].global = elementGlobal;
+ subContVol[0].inner = true;
+ subContVol[0].volume = elementVolume;
+
+ // initialize neighbors list with self:
+ numNeighbors = 1;
+ neighbors.clear();
+ neighbors.reserve(maxNE);
+ neighbors.push_back(element);
+ }
+
+ void update(const GridView& gridView, const Element& element)
+ {
+ updateInner(element);
+
+ const Geometry geometry = element.geometry();
+
+ bool onBoundary = false;
+
+ // fill neighbor information and control volume face data:
+ for (const auto& intersection : Dune::intersections(gridView, element))
+ {
+ const auto isGeometry = intersection.geometry();
+
+ // neighbor information and inner cvf data:
+ if (intersection.neighbor())
+ {
+ numNeighbors++;
+ neighbors.push_back(intersection.outside());
+
+ int scvfIdx = numNeighbors - 2;
+ SubControlVolumeFace& scvFace = subContVolFace[scvfIdx];
+
+ scvFace.i = 0;
+ scvFace.j = scvfIdx + 1;
+
+ scvFace.ipGlobal = isGeometry.center();
+ //scvFace.ipLocal = geometry.local(scvFace.ipGlobal);
+ scvFace.normal = intersection.centerUnitOuterNormal();
+ auto di = scvFace.ipGlobal;
+ di -= elementGlobal;
+ if (scvFace.normal*di < 0)
+ scvFace.normal *= -1.0;
+ Scalar volume = isGeometry.volume();
+ scvFace.normal *= volume;
+ scvFace.area = volume;
+
+ GlobalPosition distVec = elementGlobal
+ - neighbors[scvfIdx+1].geometry().center();
+ distVec /= distVec.two_norm2();
+
+ // gradients using a two-point flux approximation
+ scvFace.numFap = 2;
+ for (unsigned int fapIdx = 0; fapIdx < scvFace.numFap; fapIdx++)
+ {
+ scvFace.grad[fapIdx] = distVec;
+ scvFace.shapeValue[fapIdx] = 0.5;
+ }
+ scvFace.grad[1] *= -1.0;
+
+ scvFace.fapIndices[0] = scvFace.i;
+ scvFace.fapIndices[1] = scvFace.j;
+
+ scvFace.fIdx = intersection.indexInInside();
+ }
+
+ // boundary cvf data
+ if (intersection.boundary())
+ {
+ onBoundary = true;
+ int bfIdx = intersection.indexInInside();
+ SubControlVolumeFace& bFace = boundaryFace[bfIdx];
+
+ bFace.ipGlobal = isGeometry.center();
+ //bFace.ipLocal = geometry.local(bFace.ipGlobal);
+ bFace.normal = intersection.centerUnitOuterNormal();
+ auto di = bFace.ipGlobal;
+ di -= elementGlobal;
+ if (bFace.normal*di < 0)
+ bFace.normal *= -1.0;
+ Scalar volume = isGeometry.volume();
+ bFace.normal *= volume;
+ bFace.area = volume;
+ bFace.i = 0;
+ bFace.j = 0;
+
+ GlobalPosition distVec = elementGlobal - bFace.ipGlobal;
+ distVec /= distVec.two_norm2();
+
+ // gradients using a two-point flux approximation
+ bFace.numFap = 2;
+ for (unsigned int fapIdx = 0; fapIdx < bFace.numFap; fapIdx++)
+ {
+ bFace.grad[fapIdx] = distVec;
+ bFace.shapeValue[fapIdx] = 0.5;
+ }
+ bFace.grad[1] *= -1.0;
+
+ bFace.fapIndices[0] = bFace.i;
+ bFace.fapIndices[1] = bFace.j;
+ }
+ }
+
+ // set the number of inner-domain subcontrolvolume faces
+ numScvf = numNeighbors - 1;
+
+ // treat elements on the boundary
+ if (onBoundary)
+ {
+ for (int bfIdx = 0; bfIdx < maxBF; bfIdx++)
+ {
+ SubControlVolumeFace& bFace = boundaryFace[bfIdx];
+ bFace.j = numNeighbors + bfIdx;
+ bFace.fapIndices[1] = bFace.j;
+ neighbors.push_back(element);
+ }
+ }
+ }
+};
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh b/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7a56d40abef95467e7e2a06588cd8599457dea31
--- /dev/null
+++ b/dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh
@@ -0,0 +1,278 @@
+// -*- 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 This file contains the data which is required to calculate
+ * volume fluxes of fluid phases over a face of a finite volume by means
+ * of the Darcy approximation.
+ *
+ */
+#ifndef DUMUX_CP_DARCY_FLUX_VARIABLES_HH
+#define DUMUX_CP_DARCY_FLUX_VARIABLES_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
+NEW_PROP_TAG(SpatialParams);
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(ProblemEnableGravity);
+}
+
+/*!
+ * \ingroup ImplicitFluxVariables
+ * \brief Evaluates the normal component of the Darcy velocity
+ * on a (sub)control volume face.
+ */
+template <class TypeTag>
+class CpDarcyFluxVariables
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ enum { dim = GridView::dimension} ;
+ enum { dimWorld = GridView::dimensionworld} ;
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+ typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param fIdx The local index of the SCV (sub-control-volume) face
+ * \param elemVolVars The volume variables of the current element
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ CpDarcyFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int fIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ : fvGeometry_(fvGeometry), faceIdx_(fIdx), onBoundary_(onBoundary)
+ {
+ mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
+ calculateVolumeFlux_(problem, element, elemVolVars);
+ }
+
+public:
+ /*!
+ * \brief Return the volumetric flux over a face of a given phase.
+ *
+ * This is the calculated velocity multiplied by the unit normal
+ * and the area of the face. face().normal has already the
+ * magnitude of the area.
+ *
+ * \param phaseIdx index of the phase
+ */
+ Scalar volumeFlux(const unsigned int phaseIdx) const
+ { return volumeFlux_[phaseIdx]; }
+
+ /*!
+ * \brief Return the velocity of a given phase.
+ *
+ * This is the full velocity vector on the
+ * face (without being multiplied with normal).
+ *
+ * \param phaseIdx index of the phase
+ */
+ GlobalPosition velocity(const unsigned int phaseIdx) const
+ { return velocity_[phaseIdx] ; }
+
+ /*!
+ * \brief Return the local index of the downstream control volume
+ * for a given phase.
+ *
+ * \param phaseIdx index of the phase
+ */
+ const unsigned int downstreamIdx(const unsigned phaseIdx) const
+ { return downstreamIdx_[phaseIdx]; }
+
+ /*!
+ * \brief Return the local index of the upstream control volume
+ * for a given phase.
+ *
+ * \param phaseIdx index of the phase
+ */
+ const unsigned int upstreamIdx(const unsigned phaseIdx) const
+ { return upstreamIdx_[phaseIdx]; }
+
+ /*!
+ * \brief Return the SCV (sub-control-volume) face. This may be either
+ * a face within the element or a face on the element boundary,
+ * depending on the value of onBoundary_.
+ */
+ const SCVFace &face() const
+ {
+ if (onBoundary_)
+ return fvGeometry_.boundaryFace[faceIdx_];
+ else
+ return fvGeometry_.subContVolFace[faceIdx_];
+ }
+
+protected:
+ /*!
+ * \brief Actual calculation of the volume flux.
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemVolVars The volume variables of the current element
+ */
+ void calculateVolumeFlux_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ // calculate the transmissibilities
+ const SpatialParams &spatialParams = problem.spatialParams();
+
+ const Element& elementI = fvGeometry_.neighbors[face().i];
+ FVElementGeometry fvGeometryI;
+ fvGeometryI.subContVol[0].global = elementI.geometry().center();
+ auto ki = spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0);
+ Dune::FieldVector<Scalar, dimWorld> kin;
+ ki.mv(face().normal, kin);
+ kin /= face().area;
+ auto di = face().ipGlobal;
+ di -= elementI.geometry().center();
+ auto ti = std::abs(di*kin*face().area/(2*di.two_norm2()));
+
+ auto tij = ti;
+ if (!onBoundary_)
+ {
+ const Element& elementJ = fvGeometry_.neighbors[face().j];
+ FVElementGeometry fvGeometryJ;
+ fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+ auto kj = spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0);
+ Dune::FieldVector<Scalar, dimWorld> kjn;
+ kj.mv(face().normal, kjn);
+ kjn /= face().area;
+ auto dj = face().ipGlobal;
+ dj -= elementJ.geometry().center();
+ auto tj = std::abs(dj*kjn*face().area/(2*dj.two_norm2()));
+ tij = harmonicMean(ti, tj);
+ }
+
+ // loop over all phases
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ const auto& volVarsI = elemVolVars[face().i];
+ auto potentialI = volVarsI.pressure(phaseIdx);
+
+ const auto& volVarsJ = elemVolVars[face().j];
+ auto potentialJ = volVarsJ.pressure(phaseIdx);
+
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+ {
+ // calculate the phase density at the integration point. we
+ // only do this if the phase is present in both cells
+ Scalar SI = volVarsI.fluidState().saturation(phaseIdx);
+ Scalar SJ = volVarsJ.fluidState().saturation(phaseIdx);
+ Scalar rhoI = volVarsI.fluidState().density(phaseIdx);
+ Scalar rhoJ = volVarsJ.fluidState().density(phaseIdx);
+ Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
+ Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
+ if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
+ // doesn't matter because no wetting phase is present in
+ // both cells!
+ fI = fJ = 0.5;
+ const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
+
+ auto globalPosI = elementI.geometry().center();
+ potentialI -= density*(problem.gravityAtPos(globalPosI)*globalPosI);
+
+ if (onBoundary_)
+ {
+ potentialJ -= density*(problem.gravityAtPos(face().ipGlobal)*face().ipGlobal);
+ }
+ else
+ {
+ const Element& elementJ = fvGeometry_.neighbors[face().j];
+ auto globalPosJ = elementJ.geometry().center();
+ potentialJ -= density*(problem.gravityAtPos(globalPosJ)*globalPosJ);
+ }
+ }
+
+ auto potentialDiff = potentialI - potentialJ;
+
+ // determine the upwind direction
+ if (potentialDiff > 0)
+ {
+ upstreamIdx_[phaseIdx] = face().i;
+ downstreamIdx_[phaseIdx] = face().j;
+ }
+ else
+ {
+ upstreamIdx_[phaseIdx] = face().j;
+ downstreamIdx_[phaseIdx] = face().i;
+ }
+
+ // obtain the upwind volume variables
+ const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
+ const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
+
+ // set the volume flux
+ volumeFlux_[phaseIdx] = tij*potentialDiff;
+ volumeFlux_[phaseIdx] *= mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
+ + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx);
+
+ velocity_[phaseIdx] = face().normal;
+ velocity_[phaseIdx] *= volumeFlux_[phaseIdx]/face().area;
+ } // over loop all phases
+ }
+
+ const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
+ const unsigned int faceIdx_; //!< The index of the sub control volume face
+ const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
+ unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
+ Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
+ GlobalPosition velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
+ Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
+};
+
+} // end namespace
+
+#endif // DUMUX_CP_DARCY_FLUX_VARIABLES_HH
diff --git a/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh b/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh
index 25cfbc24f2c9ec543d268bc800f3f195c9503a58..d13d78e0014a0378038ae3da12ecb8f23ac79f95 100644
--- a/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh
+++ b/test/porousmediumflow/2p/implicit/cc2pcornerpointproblem.hh
@@ -28,9 +28,9 @@
#include <dune/grid/CpGrid.hpp>
#include <dumux/io/cpgridcreator.hh>
-#include <dumux/implicit/cornerpoint/cpelementvolumevariables.hh>
-#include <dumux/implicit/cornerpoint/cpfvelementgeometry.hh>
-#include <dumux/implicit/cornerpoint/cpdarcyfluxvariables.hh>
+#include <dumux/implicit/cornerpoint/elementvolumevariables.hh>
+#include <dumux/implicit/cornerpoint/fvelementgeometry.hh>
+#include <dumux/porousmediumflow/implicit/cpdarcyfluxvariables.hh>
#include "cc2pcornerpointspatialparams.hh"