diff --git a/dumux/decoupled/1p/1pproperties.hh b/dumux/decoupled/1p/1pproperties.hh
index 5987dab68d2d959728bda4ba36ce193bf7bf8ca2..b9552d730272024d0150c066ff584a9f6acab42e 100644
--- a/dumux/decoupled/1p/1pproperties.hh
+++ b/dumux/decoupled/1p/1pproperties.hh
@@ -30,6 +30,8 @@
#ifndef DUMUX_1PPROPERTIES_HH
#define DUMUX_1PPROPERTIES_HH
+#define OnePModel
+
//Dune-includes
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
diff --git a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
index 307542b226ce534179bbd6912fa44bbd4be2f135..af492500caca34e2922d4517e2611410654af13d 100644
--- a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
+++ b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
@@ -40,13 +40,15 @@ namespace Dumux
* @tparam TypeTag The Type Tag
* @tparam Implementation The Problem implementation
*/
-template<class TypeTag, class Implementation>
+template<class TypeTag>
class DiffusionProblem1P: public OneModelProblem<TypeTag>
{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Implementation;
typedef OneModelProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GridView::Grid Grid;typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
@@ -54,6 +56,8 @@ class DiffusionProblem1P: public OneModelProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+
enum
{
dim = Grid::dimension, dimWorld = Grid::dimensionworld
@@ -68,41 +72,8 @@ public:
* \param gridView The grid view
* \param verbose Output flag for the time manager.
*/
- DiffusionProblem1P(const GridView &gridView, bool verbose = true)
- DUNE_DEPRECATED // use DiffusionProblem1P(TimeManager&, const GridView&)
- : ParentType(gridView, verbose), gravity_(0)
- {
- spatialParameters_ = new SpatialParameters(gridView);
- newSpatialParams_ = true;
- gravity_ = 0;
- if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
- gravity_[dim - 1] = -9.81;
- }
- /*!
- * \brief The constructor
- *
- * \param gridView The grid view
- * \param spatialParameters SpatialParameters instantiation
- * \param verbose Output flag for the time manager.
- */
- DiffusionProblem1P(const GridView &gridView, SpatialParameters &spatialParameters, bool verbose = true)
- DUNE_DEPRECATED // use DiffusionProblem1P(TimeManager&, const GridView&, SpatialParameters&)
- : ParentType(gridView, verbose), gravity_(0), spatialParameters_(&spatialParameters)
- {
- newSpatialParams_ = false;
- gravity_ = 0;
- if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
- gravity_[dim - 1] = -9.81;
- }
-
- /*!
- * \brief The constructor
- *
- * \param gridView The grid view
- * \param verbose Output flag for the time manager.
- */
- DiffusionProblem1P(TimeManager &timeManager, const GridView &gridView) :
- ParentType(timeManager, gridView), gravity_(0)
+ DiffusionProblem1P(const GridView &gridView)
+ : ParentType(gridView), gravity_(0)
{
spatialParameters_ = new SpatialParameters(gridView);
newSpatialParams_ = true;
@@ -117,8 +88,8 @@ public:
* \param spatialParameters SpatialParameters instantiation
* \param verbose Output flag for the time manager.
*/
- DiffusionProblem1P(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters) :
- ParentType(timeManager, gridView), gravity_(0), spatialParameters_(&spatialParameters)
+ DiffusionProblem1P(const GridView &gridView, SpatialParameters &spatialParameters)
+ : ParentType(gridView), gravity_(0), spatialParameters_(&spatialParameters)
{
newSpatialParams_ = false;
gravity_ = 0;
@@ -159,13 +130,52 @@ public:
/*!
* \brief Returns the temperature within the domain.
*
- * This method MUST be overwritten by the actual problem.
+ * \param element The element
+ *
+ */
+ Scalar temperature(const Element& element) const
+ {
+ return asImp_().temperatureAtPos(element.geometry().center());
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * \param globalPos The position of the center of an element
+ *
*/
- Scalar temperature() const
+ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
- return this->asImp_()->temperature();
+ // Throw an exception (there is no initial condition)
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide "
+ "a temperatureAtPos() method.");
+ }
+
+ /*!
+ * \brief Returns the reference pressure for evaluation of constitutive relations.
+ *
+ * \param element The element
+ *
+ */
+ Scalar referencePressure(const Element& element) const
+ {
+ return asImp_().referencePressureAtPos(element.geometry().center());
+ }
+
+ /*!
+ * \brief Returns the reference pressure for evaluation of constitutive relations.
+ *
+ * \param globalPos The position of the center of an element
+ *
+ */
+ Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
+ {
+ // Throw an exception (there is no initial condition)
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide "
+ "a referencePressureAtPos() method.");
}
- ;
/*!
* \brief Returns the acceleration due to gravity.
@@ -197,6 +207,15 @@ public:
// \}
private:
+private:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this); }
+
+ //! \copydoc Dumux::IMPETProblem::asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this); }
+
GlobalPosition gravity_;
// fluids and material properties
diff --git a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
index 67510305e85da6e4b3e09fe98500c537deb2c531..35245c1f12f89813ebe600402e7958a71b3b8480 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
@@ -295,9 +295,6 @@ void FVPressure1P<TypeTag>::assemble(bool first)
f_[globalIdxI] = volume * source;
- // get absolute permeability
- FieldMatrix permeabilityI(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
-
int isIndex = 0;
IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt);
@@ -326,27 +323,14 @@ void FVPressure1P<TypeTag>::assemble(bool first)
// compute distance between cell centers
Scalar dist = distVec.two_norm();
- FieldMatrix permeabilityJ = problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor,
- *neighborPointer);
-
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
- // harmonic mean of permeability
- for (int x = 0; x < dim; x++)
- {
- meanPermeability[x][x] = 2 * permeabilityI[x][x] * permeabilityJ[x][x] / (permeabilityI[x][x]
- + permeabilityJ[x][x]);
- for (int y = 0; y < dim; y++)
- {
- if (x != y)
- {//use arithmetic mean for the off-diagonal entries to keep the tensor property!
- meanPermeability[x][y] = 0.5 * (permeabilityI[x][y] + permeabilityJ[x][y]);
- }
- }
- }
+ problem_.spatialParameters().meanK(meanPermeability,
+ problem_.spatialParameters().intrinsicPermeability(*eIt),
+ problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
- Dune::FieldVector < Scalar, dim > permeability(0);
+ Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
permeability/=viscosityI;
@@ -413,17 +397,24 @@ void FVPressure1P<TypeTag>::assemble(bool first)
// center of face in global coordinates
const GlobalPosition& globalPosFace = isIt->geometry().center();
- Dune::FieldVector < Scalar, dimWorld > distVec(globalPosFace - globalPos);
- Scalar dist = distVec.two_norm();
-
//get boundary condition for boundary face center
BoundaryConditions::Flags bctype = problem_.bctype(globalPosFace, *isIt);
if (bctype == BoundaryConditions::dirichlet)
{
+ Dune::FieldVector < Scalar, dimWorld > distVec(globalPosFace - globalPos);
+ Scalar dist = distVec.two_norm();
+
+ //permeability vector at boundary
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+
+ problem_.spatialParameters().meanK(meanPermeability,
+ problem_.spatialParameters().intrinsicPermeability(*eIt));
+
//permeability vector at boundary
Dune::FieldVector < Scalar, dim > permeability(0);
- permeabilityI.mv(unitOuterNormal, permeability);
+ meanPermeability.mv(unitOuterNormal, permeability);
permeability/= viscosityI;
diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index ec792d6832ac5126a3da226a0ba5ebffdff8b2f0..3276a7d3d2b131bc0341f397ea5871b9028482e9 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -189,9 +189,6 @@ void FVVelocity1P<TypeTag>::calculateVelocity()
Dune::FieldVector<Scalar,dimWorld> unitOuterNormal = isIt->centerUnitOuterNormal();
- // get absolute permeability
- FieldMatrix permeabilityI(this->problem().spatialParameters().intrinsicPermeability(globalPos, *eIt));
-
// handle interior face
if (isIt->neighbor())
{
@@ -212,30 +209,18 @@ void FVVelocity1P<TypeTag>::calculateVelocity()
// compute distance between cell centers
Scalar dist = distVec.two_norm();
- // get absolute permeability
- FieldMatrix permeabilityJ(this->problem().spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
-
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
-// // harmonic mean of permeability
- for (int x = 0; x < dim; x++)
- {
- meanPermeability[x][x] = 2 * permeabilityI[x][x] * permeabilityJ[x][x] / (permeabilityI[x][x]
- + permeabilityJ[x][x]);
- for (int y = 0; y < dim; y++)
- {
- if (x != y)
- {//use arithmetic mean for the off-diagonal entries to keep the tensor property!
- meanPermeability[x][y] = 0.5 * (permeabilityI[x][y] + permeabilityJ[x][y]);
- }
- }
- }
-
- Dune::FieldVector<Scalar,dim> permeability(0);
- meanPermeability.mv(unitOuterNormal,permeability);
+ this->problem().spatialParameters().meanK(meanPermeability,
+ this->problem().spatialParameters().intrinsicPermeability(*eIt),
+ this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
permeability/=viscosityI;
+
// std::cout<<"Mean Permeability / Viscosity: "<<meanPermeability<<std::endl;
Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
@@ -284,22 +269,29 @@ void FVVelocity1P<TypeTag>::calculateVelocity()
//get boundary type
BoundaryConditions::Flags bcType = this->problem().bctype(globalPosFace, *isIt);
- // cell center in global coordinates
- GlobalPosition globalPos = eIt->geometry().center();
+ if (bcType == BoundaryConditions::dirichlet)
+ {
+ // cell center in global coordinates
+ GlobalPosition globalPos = eIt->geometry().center();
- // distance vector between barycenters
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosFace - globalPos;
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar,dimWorld> distVec = globalPosFace - globalPos;
- // compute distance between cell centers
- Scalar dist = distVec.two_norm();
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
- //multiply with normal vector at the boundary
- Dune::FieldVector<Scalar,dim> permeability(0);
- permeabilityI.mv(unitOuterNormal, permeability);
- permeability/=viscosityI;
+ //permeability vector at boundary
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+
+ this->problem().spatialParameters().meanK(meanPermeability,
+ this->problem().spatialParameters().intrinsicPermeability(*eIt));
+
+ //multiply with normal vector at the boundary
+ Dune::FieldVector<Scalar,dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
+ permeability/=viscosityI;
- if (bcType == BoundaryConditions::dirichlet)
- {
Scalar pressBound = this->problem().dirichlet(globalPosFace, *isIt);
//calculate the gravity term
diff --git a/dumux/decoupled/common/onemodelproblem.hh b/dumux/decoupled/common/onemodelproblem.hh
index a6e19bd2bec2b54993c93b95efb3679124271ec5..1e8d044f3f66e34a102c9de1c31d7f278162a15c 100644
--- a/dumux/decoupled/common/onemodelproblem.hh
+++ b/dumux/decoupled/common/onemodelproblem.hh
@@ -92,7 +92,6 @@ public:
* \tparam verbose Output level for Dumux::TimeManager
*/
OneModelProblem(const GridView &gridView, bool verbose = true)
- DUNE_DEPRECATED // use OneModelProblem(TimeManager&, const GridView&)
: gridView_(gridView),
bboxMin_(std::numeric_limits<double>::max()),
bboxMax_(-std::numeric_limits<double>::max()),