diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 434f45874496dd82257997dfc64edd991fcdf9d5..185b549560675b0dbb4ca9a491e19f2a0bf32685 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/dumux/boxmodels/2p/2pmodel.hh
@@ -74,7 +74,7 @@ namespace Dumux
* default, the model uses \f$p_w\f$ and \f$S_n\f$.
*/
template<class TypeTag >
-class TwoPModel : public BoxModel<TypeTag>
+class TwoPModel : public GET_PROP_TYPE(TypeTag, BaseModel)
{
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
@@ -133,22 +133,29 @@ public:
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
- // create the required scalar fields
+ // get the number of degrees of freedom and the number of vertices
+ unsigned numDofs = this->numDofs();
unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField *pW = writer.allocateManagedBuffer(numVertices);
- ScalarField *pN = writer.allocateManagedBuffer(numVertices);
- ScalarField *pC = writer.allocateManagedBuffer(numVertices);
- ScalarField *Sw = writer.allocateManagedBuffer(numVertices);
- ScalarField *Sn = writer.allocateManagedBuffer(numVertices);
- ScalarField *rhoW = writer.allocateManagedBuffer(numVertices);
- ScalarField *rhoN = writer.allocateManagedBuffer(numVertices);
- ScalarField *mobW = writer.allocateManagedBuffer(numVertices);
- ScalarField *mobN = writer.allocateManagedBuffer(numVertices);
- ScalarField *poro = writer.allocateManagedBuffer(numVertices);
- ScalarField *Te = writer.allocateManagedBuffer(numVertices);
- ScalarField *cellNum =writer.allocateManagedBuffer (numVertices);
- VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numVertices);
- VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numVertices);
+
+ // velocity output currently only works for vertex data
+ if (numDofs != numVertices)
+ velocityOutput = false;
+
+ // create the required scalar fields
+ ScalarField *pW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pC = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sw = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sn = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Te = writer.allocateManagedBuffer(numDofs);
+ ScalarField *cellNum =writer.allocateManagedBuffer (numDofs);
+ VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
if(velocityOutput) // check if velocity output is demanded
{
@@ -181,172 +188,210 @@ public:
fvElemGeom.update(this->gridView_(), *elemIt);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ if (numDofs == numElements) // element data
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
- volVars.update(sol[globalIdx],
+ volVars.update(sol[idx],
this->problem_(),
*elemIt,
fvElemGeom,
- i,
+ /*scvIdx=*/0,
false);
- (*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
- (*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
- (*pC)[globalIdx] = volVars.capillaryPressure();
- (*Sw)[globalIdx] = volVars.saturation(wPhaseIdx);
- (*Sn)[globalIdx] = volVars.saturation(nPhaseIdx);
- (*rhoW)[globalIdx] = volVars.density(wPhaseIdx);
- (*rhoN)[globalIdx] = volVars.density(nPhaseIdx);
- (*mobW)[globalIdx] = volVars.mobility(wPhaseIdx);
- (*mobN)[globalIdx] = volVars.mobility(nPhaseIdx);
- (*poro)[globalIdx] = volVars.porosity();
- (*Te)[globalIdx] = volVars.temperature();
- if(velocityOutput)
- {
- (*cellNum)[globalIdx] += 1;
- }
- };
-
- if(velocityOutput)
+ (*pW)[idx] = volVars.pressure(wPhaseIdx);
+ (*pN)[idx] = volVars.pressure(nPhaseIdx);
+ (*pC)[idx] = volVars.capillaryPressure();
+ (*Sw)[idx] = volVars.saturation(wPhaseIdx);
+ (*Sn)[idx] = volVars.saturation(nPhaseIdx);
+ (*rhoW)[idx] = volVars.density(wPhaseIdx);
+ (*rhoN)[idx] = volVars.density(nPhaseIdx);
+ (*mobW)[idx] = volVars.mobility(wPhaseIdx);
+ (*mobN)[idx] = volVars.mobility(nPhaseIdx);
+ (*poro)[idx] = volVars.porosity();
+ (*Te)[idx] = volVars.temperature();
+ }
+ else // vertex data
{
- // calculate vertex velocities
- GlobalPosition tmpVelocity[numPhases];
-
- for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ int numVerts = elemIt->template count<dim> ();
+ for (int i = 0; i < numVerts; ++i)
{
- tmpVelocity[phaseIdx] = Scalar(0.0);
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
+
+ (*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
+ (*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
+ (*pC)[globalIdx] = volVars.capillaryPressure();
+ (*Sw)[globalIdx] = volVars.saturation(wPhaseIdx);
+ (*Sn)[globalIdx] = volVars.saturation(nPhaseIdx);
+ (*rhoW)[globalIdx] = volVars.density(wPhaseIdx);
+ (*rhoN)[globalIdx] = volVars.density(nPhaseIdx);
+ (*mobW)[globalIdx] = volVars.mobility(wPhaseIdx);
+ (*mobN)[globalIdx] = volVars.mobility(nPhaseIdx);
+ (*poro)[globalIdx] = volVars.porosity();
+ (*Te)[globalIdx] = volVars.temperature();
+ if(velocityOutput)
+ {
+ (*cellNum)[globalIdx] += 1;
+ }
}
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
- SCVVelocities scvVelocityW(8), scvVelocityN(8);
+ if(velocityOutput)
+ {
+ // calculate vertex velocities
+ GlobalPosition tmpVelocity[numPhases];
- scvVelocityW = 0;
- scvVelocityN = 0;
+ for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ tmpVelocity[phaseIdx] = Scalar(0.0);
+ }
- ElementVolumeVariables elemVolVars;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
+ SCVVelocities scvVelocityW(8), scvVelocityN(8);
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvElemGeom,
- false /* oldSol? */);
+ scvVelocityW = 0;
+ scvVelocityN = 0;
- for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
- {
+ ElementVolumeVariables elemVolVars;
- FluxVariables fluxDat(this->problem_(),
- *elemIt,
- fvElemGeom,
- faceIdx,
- elemVolVars);
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ false /* oldSol? */);
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
{
- // data attached to upstream and the downstream vertices
- // of the current phase
- const VolumeVariables up =
- elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
- const VolumeVariables dn =
- elemVolVars[fluxDat.downstreamIdx(phaseIdx)];
-
- // calculate the flux in the normal direction of the
- // current sub control volume face
- GlobalPosition tmpVec(0);
- fluxDat.intrinsicPermeability().mv(fluxDat.potentialGrad(phaseIdx),
- tmpVec);
- const GlobalPosition globalNormal = fluxDat.face().normal;
- const Scalar normalFlux = - (tmpVec*globalNormal);
-
- // local position of integration point
- const Dune::FieldVector<Scalar, dim>& localPosIP = fvElemGeom.subContVolFace[faceIdx].ipLocal;
-
- // Transformation of the global normal vector to normal vector in the reference element
- const Dune::FieldMatrix<CoordScalar, dim, dim> jacobianT1 = elemIt->geometry().jacobianTransposed(localPosIP);
-
- GlobalPosition localNormal(0);
- jacobianT1.mv(globalNormal, localNormal);
- // note only works for cubes
- const Scalar localArea = pow(2,-(dim-1));
-
- localNormal /= localNormal.two_norm();
-
- // Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolumeface
- // in the reference element.
- const Scalar massUpwindWeight = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
- PhasesVector q;
- q[phaseIdx] = normalFlux
+ FluxVariables fluxDat(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ faceIdx,
+ elemVolVars);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables up =
+ elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
+ const VolumeVariables dn =
+ elemVolVars[fluxDat.downstreamIdx(phaseIdx)];
+
+ // calculate the flux in the normal direction of the
+ // current sub control volume face
+ GlobalPosition tmpVec(0);
+ fluxDat.intrinsicPermeability().mv(fluxDat.potentialGrad(phaseIdx),
+ tmpVec);
+ const GlobalPosition globalNormal = fluxDat.face().normal;
+ const Scalar normalFlux = - (tmpVec*globalNormal);
+
+ // local position of integration point
+ const Dune::FieldVector<Scalar, dim>& localPosIP = fvElemGeom.subContVolFace[faceIdx].ipLocal;
+
+ // Transformation of the global normal vector to normal vector in the reference element
+ const Dune::FieldMatrix<CoordScalar, dim, dim> jacobianT1 = elemIt->geometry().jacobianTransposed(localPosIP);
+
+ GlobalPosition localNormal(0);
+ jacobianT1.mv(globalNormal, localNormal);
+ // note only works for cubes
+ const Scalar localArea = pow(2,-(dim-1));
+
+ localNormal /= localNormal.two_norm();
+
+ // Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolumeface
+ // in the reference element.
+ const Scalar massUpwindWeight = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
+ PhasesVector q;
+ q[phaseIdx] = normalFlux
* (massUpwindWeight
* up.mobility(phaseIdx)
+ (1- massUpwindWeight)
* dn.mobility(phaseIdx)) / localArea;
- // transform the normal Darcy velocity into a vector
- tmpVelocity[phaseIdx] = localNormal;
- tmpVelocity[phaseIdx] *= q[phaseIdx];
-
- if(phaseIdx == wPhaseIdx){
- scvVelocityW[fluxDat.face().i] += tmpVelocity[phaseIdx];
- scvVelocityW[fluxDat.face().j] += tmpVelocity[phaseIdx];
- }
- else if(phaseIdx == nPhaseIdx){
- scvVelocityN[fluxDat.face().i] += tmpVelocity[phaseIdx];
- scvVelocityN[fluxDat.face().j] += tmpVelocity[phaseIdx];
- }
- }
- }
- typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
- const Dune::FieldVector<Scalar, dim> &localPos
- = ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
-
- // get the transposed Jacobian of the element mapping
- const Dune::FieldMatrix<CoordScalar, dim, dim> &jacobianT2
- = elemIt->geometry().jacobianTransposed(localPos);
-
- // transform vertex velocities from local to global coordinates
- for (int i = 0; i < numVerts; ++i)
- {
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
- // calculate the subcontrolvolume velocity by the Piola transformation
- Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
-
- jacobianT2.mtv(scvVelocityW[i], scvVelocity);
- scvVelocity /= elemIt->geometry().integrationElement(localPos);
- // add up the wetting phase subcontrolvolume velocities for each vertex
- (*velocityW)[globalIdx] += scvVelocity;
-
- jacobianT2.mtv(scvVelocityN[i], scvVelocity);
- scvVelocity /= elemIt->geometry().integrationElement(localPos);
- // add up the nonwetting phase subcontrolvolume velocities for each vertex
- (*velocityN)[globalIdx] += scvVelocity;
+ // transform the normal Darcy velocity into a vector
+ tmpVelocity[phaseIdx] = localNormal;
+ tmpVelocity[phaseIdx] *= q[phaseIdx];
+
+ if(phaseIdx == wPhaseIdx){
+ scvVelocityW[fluxDat.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityW[fluxDat.face().j] += tmpVelocity[phaseIdx];
+ }
+ else if(phaseIdx == nPhaseIdx){
+ scvVelocityN[fluxDat.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityN[fluxDat.face().j] += tmpVelocity[phaseIdx];
+ }
+ }
+ }
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::FieldVector<Scalar, dim> &localPos
+ = ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
+
+ // get the transposed Jacobian of the element mapping
+ const Dune::FieldMatrix<CoordScalar, dim, dim> &jacobianT2
+ = elemIt->geometry().jacobianTransposed(localPos);
+
+ // transform vertex velocities from local to global coordinates
+ for (int i = 0; i < numVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ // calculate the subcontrolvolume velocity by the Piola transformation
+ Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
+
+ jacobianT2.mtv(scvVelocityW[i], scvVelocity);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
+ // add up the wetting phase subcontrolvolume velocities for each vertex
+ (*velocityW)[globalIdx] += scvVelocity;
+
+ jacobianT2.mtv(scvVelocityN[i], scvVelocity);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
+ // add up the nonwetting phase subcontrolvolume velocities for each vertex
+ (*velocityN)[globalIdx] += scvVelocity;
+ }
}
}
}
- if(velocityOutput)
+
+ if (numDofs == numElements) // element data
+ {
+ writer.attachCellData(*Sn, "Sn");
+ writer.attachCellData(*Sw, "Sw");
+ writer.attachCellData(*pN, "pn");
+ writer.attachCellData(*pW, "pw");
+ writer.attachCellData(*pC, "pc");
+ writer.attachCellData(*rhoW, "rhoW");
+ writer.attachCellData(*rhoN, "rhoN");
+ writer.attachCellData(*mobW, "mobW");
+ writer.attachCellData(*mobN, "mobN");
+ writer.attachCellData(*poro, "porosity");
+ writer.attachCellData(*Te, "temperature");
+ }
+ else // vertex data
+ {
+ writer.attachVertexData(*Sn, "Sn");
+ writer.attachVertexData(*Sw, "Sw");
+ writer.attachVertexData(*pN, "pn");
+ writer.attachVertexData(*pW, "pw");
+ writer.attachVertexData(*pC, "pc");
+ writer.attachVertexData(*rhoW, "rhoW");
+ writer.attachVertexData(*rhoN, "rhoN");
+ writer.attachVertexData(*mobW, "mobW");
+ writer.attachVertexData(*mobN, "mobN");
+ writer.attachVertexData(*poro, "porosity");
+ writer.attachVertexData(*Te, "temperature");
+ if(velocityOutput) // check if velocity output is demanded
{
// divide the vertex velocities by the number of adjacent scvs i.e. cells
- for(int globalIdx = 0; globalIdx<numVertices; ++globalIdx){
- (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
- (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
+ for(int globalIdx = 0; globalIdx < numVertices; ++globalIdx){
+ (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
+ (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
}
+ writer.attachVertexData(*velocityW, "velocityW", dim);
+ writer.attachVertexData(*velocityN, "velocityN", dim);
}
-
- writer.attachVertexData(*Sn, "Sn");
- writer.attachVertexData(*Sw, "Sw");
- writer.attachVertexData(*pN, "pn");
- writer.attachVertexData(*pW, "pw");
- writer.attachVertexData(*pC, "pc");
- writer.attachVertexData(*rhoW, "rhoW");
- writer.attachVertexData(*rhoN, "rhoN");
- writer.attachVertexData(*mobW, "mobW");
- writer.attachVertexData(*mobN, "mobN");
- writer.attachVertexData(*poro, "porosity");
- writer.attachVertexData(*Te, "temperature");
- if(velocityOutput) // check if velocity output is demanded
- {
- writer.attachVertexData(*velocityW, "velocityW", dim);
- writer.attachVertexData(*velocityN, "velocityN", dim);
}
writer.attachCellData(*rank, "process rank");
}
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index ca54b0c383bca0680d9c5f2dbc47c4254eda81f1..e5a24b9aab8aae9fe9ff426e75a0fd20cebbafab 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -64,7 +64,8 @@ NEW_PROP_TAG(GridView); //!< The type of the grid view
NEW_PROP_TAG(FVElementGeometry); //! The type of the finite-volume geometry in the box scheme
NEW_PROP_TAG(Problem); //!< The type of the problem
-NEW_PROP_TAG(Model); //!< The type of the discretization
+NEW_PROP_TAG(BaseModel); //!< The type of the base class of the model
+NEW_PROP_TAG(Model); //!< The type of the model
NEW_PROP_TAG(NumEq); //!< Number of equations in the system of PDEs
NEW_PROP_TAG(BaseLocalResidual); //!< The type of the base class of the local residual
NEW_PROP_TAG(LocalResidual); //!< The type of the local residual function
diff --git a/dumux/boxmodels/common/boxpropertydefaults.hh b/dumux/boxmodels/common/boxpropertydefaults.hh
index 6fab788fa88e489957a7271b86e27e34f2172737..2031c43ec6de585d19a7167462ceb442bc3c2f62 100644
--- a/dumux/boxmodels/common/boxpropertydefaults.hh
+++ b/dumux/boxmodels/common/boxpropertydefaults.hh
@@ -35,6 +35,7 @@
#include <dumux/nonlinear/newtoncontroller.hh>
#include "boxassembler.hh"
+#include "boxmodel.hh"
#include "boxfvelementgeometry.hh"
#include "boxelementboundarytypes.hh"
#include "boxlocaljacobian.hh"
@@ -51,6 +52,10 @@
namespace Dumux {
+// forward declaration
+template<class TypeTag>
+class BoxModel;
+
namespace Properties {
//////////////////////////////////////////////////////////////////
// Some defaults for very fundamental properties
@@ -98,6 +103,9 @@ SET_TYPE_PROP(BoxModel, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper)
//! Set the BaseLocalResidual to BoxLocalResidual
SET_TYPE_PROP(BoxModel, BaseLocalResidual, Dumux::BoxLocalResidual<TypeTag>);
+//! Set the BaseModel to BoxModel
+SET_TYPE_PROP(BoxModel, BaseModel, Dumux::BoxModel<TypeTag>);
+
//! The local jacobian operator for the box scheme
SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);