From bb0e3a3077d9ff8bb32121fd8698f5efe6e1c0dc Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Mon, 10 Dec 2012 14:51:36 +0000
Subject: [PATCH] implicit branch: add attachDofData method to the
VTKMultiWriter to uniformly treat vertex and cell data. Simplify
addOutputVtkFields methods of the models correspondingly. Get the global
indices from the dofMapper instead of vertexMapper or elementMapper whereever
possible, avoiding corresponding if-statements.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9792 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/1p/1pmodel.hh | 22 +--
dumux/implicit/1p2c/1p2cmodel.hh | 71 +++-----
dumux/implicit/2p/2pmodel.hh | 80 ++++-----
dumux/implicit/2p2c/2p2cmodel.hh | 153 +++++++-----------
dumux/implicit/2p2c/2p2cvolumevariables.hh | 10 +-
dumux/implicit/3p3c/3p3cmodel.hh | 110 ++++---------
dumux/implicit/3p3c/3p3cvolumevariables.hh | 7 +-
dumux/implicit/co2/co2model.hh | 28 ++--
dumux/implicit/co2/co2volumevariables.hh | 13 +-
dumux/implicit/common/implicitmodel.hh | 15 +-
.../mpnc/energy/mpncvtkwriterenergy.hh | 20 +--
dumux/implicit/mpnc/mass/mpncvtkwritermass.hh | 12 +-
dumux/implicit/mpnc/mpncvtkwritercommon.hh | 7 +-
dumux/implicit/richards/richardsmodel.hh | 46 ++----
.../richards/richardsnewtoncontroller.hh | 7 +-
dumux/io/vtkmultiwriter.hh | 16 ++
16 files changed, 223 insertions(+), 394 deletions(-)
diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index 40b419d33f..01603e78c2 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -60,6 +60,7 @@ class OnePBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
enum { dim = GridView::dimension };
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -98,11 +99,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox)
- globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*elemIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
@@ -113,23 +110,14 @@ public:
const SpatialParams &spatialParams = this->problem_().spatialParams();
(*p)[globalIdx] = volVars.pressure();
- (*K)[globalIdx]= spatialParams.intrinsicPermeability(*elemIt,
+ (*K)[globalIdx] = spatialParams.intrinsicPermeability(*elemIt,
fvGeometry,
scvIdx);
}
}
- if (isBox)
- {
- writer.attachVertexData(*p, "p");
- writer.attachVertexData(*K, "K");
- }
- else
- {
- writer.attachCellData(*p, "p");
- writer.attachCellData(*K, "K");
- }
-
+ writer.attachDofData(*p, "p", isBox);
+ writer.attachDofData(*K, "K", isBox);
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index d565a1fa1a..75d7264638 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -91,6 +91,7 @@ class OnePTwoCBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, dim> DimVector;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -175,11 +176,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox)
- globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*elemIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
@@ -245,7 +242,7 @@ public:
//use vertical faces (horizontal noraml vector) to calculate vx
//in case of heterogeneities it seams to be better to define intrinisc permeability elementwise
- if(xDir > yDir)//(fluxVars.face().normal[0] > 1e-10 || fluxVars.face().normal[0] < -1e-10)// (xDir > yDir)
+ if (xDir > yDir)//(fluxVars.face().normal[0] > 1e-10 || fluxVars.face().normal[0] < -1e-10)// (xDir > yDir)
{
// get darcy velocity
//calculate (v n) n/A
@@ -309,50 +306,28 @@ public:
}
}
- if (isBox)
- {
- writer.attachVertexData(pressure, "P");
- writer.attachVertexData(delp, "delp");
- if (velocityOutput)
- {
- writer.attachVertexData(velocityX, "Vx");
- writer.attachVertexData(velocityY, "Vy");
- if (dim > 2)
- writer.attachVertexData(velocityZ, "Vz");
- }
- char nameMoleFraction0[42], nameMoleFraction1[42];
- snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
- snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
- writer.attachVertexData(moleFraction0, nameMoleFraction0);
- writer.attachVertexData(moleFraction1, nameMoleFraction1);
-
- char nameMassFraction0[42], nameMassFraction1[42];
- snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
- snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
- writer.attachVertexData(massFraction0, nameMassFraction0);
- writer.attachVertexData(massFraction1, nameMassFraction1);
- writer.attachVertexData(rho, "rho");
- writer.attachVertexData(mu, "mu");
- }
- else
+ writer.attachDofData(pressure, "P", isBox);
+ writer.attachDofData(delp, "delp", isBox);
+ if (velocityOutput)
{
- writer.attachCellData(pressure, "P");
- writer.attachCellData(delp, "delp");
- char nameMoleFraction0[42], nameMoleFraction1[42];
- snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
- snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
- writer.attachCellData(moleFraction0, nameMoleFraction0);
- writer.attachCellData(moleFraction1, nameMoleFraction1);
-
- char nameMassFraction0[42], nameMassFraction1[42];
- snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
- snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
- writer.attachCellData(massFraction0, nameMassFraction0);
- writer.attachCellData(massFraction1, nameMassFraction1);
- writer.attachCellData(rho, "rho");
- writer.attachCellData(mu, "mu");
+ writer.attachDofData(velocityX, "Vx", isBox);
+ writer.attachDofData(velocityY, "Vy", isBox);
+ if (dim > 2)
+ writer.attachDofData(velocityZ, "Vz", isBox);
}
-
+ char nameMoleFraction0[42], nameMoleFraction1[42];
+ snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
+ snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
+ writer.attachDofData(moleFraction0, nameMoleFraction0, isBox);
+ writer.attachDofData(moleFraction1, nameMoleFraction1, isBox);
+
+ char nameMassFraction0[42], nameMassFraction1[42];
+ snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
+ snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
+ writer.attachDofData(massFraction0, nameMassFraction0, isBox);
+ writer.attachDofData(massFraction1, nameMassFraction1, isBox);
+ writer.attachDofData(rho, "rho", isBox);
+ writer.attachDofData(mu, "mu", isBox);
writer.attachCellData(rank, "process rank");
}
diff --git a/dumux/implicit/2p/2pmodel.hh b/dumux/implicit/2p/2pmodel.hh
index 4092214183..a045ae796e 100644
--- a/dumux/implicit/2p/2pmodel.hh
+++ b/dumux/implicit/2p/2pmodel.hh
@@ -99,6 +99,7 @@ class TwoPModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -140,7 +141,7 @@ public:
VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
- if(velocityOutput) // check if velocity output is demanded
+ if (velocityOutput) // check if velocity output is demanded
{
// initialize velocity fields
for (unsigned int i = 0; i < numDofs; ++i)
@@ -161,7 +162,7 @@ public:
ElementIterator elemEndIt = this->gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt)
{
- if(velocityOutput && !elemIt->geometry().type().isCube()){
+ if (velocityOutput && !elemIt->geometry().type().isCube()){
DUNE_THROW(Dune::InvalidStateException,
"Currently, velocity output only works for cubes. "
"Please set the EnableVelocityOutput property to false!");
@@ -173,11 +174,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox)
- globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*elemIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
@@ -197,13 +194,13 @@ public:
(*mobN)[globalIdx] = volVars.mobility(nPhaseIdx);
(*poro)[globalIdx] = volVars.porosity();
(*Te)[globalIdx] = volVars.temperature();
- if(velocityOutput)
+ if (velocityOutput)
{
(*cellNum)[globalIdx] += 1;
}
}
- if(velocityOutput)
+ if (velocityOutput)
{
// calculate vertex velocities
GlobalPosition tmpVelocity[numPhases];
@@ -261,11 +258,11 @@ public:
tmpVelocity[phaseIdx] = localNormal;
tmpVelocity[phaseIdx] *= q[phaseIdx];
- if(phaseIdx == wPhaseIdx){
+ if (phaseIdx == wPhaseIdx){
scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
- else if(phaseIdx == nPhaseIdx){
+ else if (phaseIdx == nPhaseIdx){
scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
@@ -280,18 +277,18 @@ public:
= elemIt->geometry().jacobianTransposed(localPos);
// transform vertex velocities from local to global coordinates
- for (int i = 0; i < fvGeometry.numSCV; ++i)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
- jacobianT2.mtv(scvVelocityW[i], scvVelocity);
+ jacobianT2.mtv(scvVelocityW[scvIdx], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the wetting phase subcontrolvolume velocities for each vertex
(*velocityW)[globalIdx] += scvVelocity;
- jacobianT2.mtv(scvVelocityN[i], scvVelocity);
+ jacobianT2.mtv(scvVelocityN[scvIdx], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the nonwetting phase subcontrolvolume velocities for each vertex
(*velocityN)[globalIdx] += scvVelocity;
@@ -299,43 +296,26 @@ public:
}
}
- if (isBox) // vertex data
+ writer.attachDofData(*Sn, "Sn", isBox);
+ writer.attachDofData(*Sw, "Sw", isBox);
+ writer.attachDofData(*pN, "pn", isBox);
+ writer.attachDofData(*pW, "pw", isBox);
+ writer.attachDofData(*pC, "pc", isBox);
+ writer.attachDofData(*rhoW, "rhoW", isBox);
+ writer.attachDofData(*rhoN, "rhoN", isBox);
+ writer.attachDofData(*mobW, "mobW", isBox);
+ writer.attachDofData(*mobN, "mobN", isBox);
+ writer.attachDofData(*poro, "porosity", isBox);
+ writer.attachDofData(*Te, "temperature", isBox);
+ if (velocityOutput) // check if velocity output is demanded
{
- 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(unsigned int globalIdx = 0; globalIdx < velocityW->size(); ++globalIdx){
- (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
- (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
- }
- writer.attachVertexData(*velocityW, "velocityW", dim);
- writer.attachVertexData(*velocityN, "velocityN", dim);
+ // divide the vertex velocities by the number of adjacent scvs i.e. cells
+ for(unsigned int globalIdx = 0; globalIdx < velocityW->size(); ++globalIdx){
+ (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
+ (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
}
- }
- else // cell 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");
+ writer.attachDofData(*velocityW, "velocityW", isBox, dim);
+ writer.attachDofData(*velocityN, "velocityN", isBox, dim);
}
writer.attachCellData(*rank, "process rank");
}
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index 19e9475954..d7c295649d 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -135,6 +135,7 @@ class TwoPTwoCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -154,11 +155,11 @@ public:
// check, if velocity output can be used (works only for cubes so far)
velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
- ElementIterator eIt = this->gridView_().template begin<0>();
- ElementIterator eEndIt = this->gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt)
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt)
{
- if(eIt->geometry().type().isCube() == false){
+ if (elemIt->geometry().type().isCube() == false){
velocityOutput_ = false;
}
@@ -166,8 +167,8 @@ public:
{
velocityOutput_ = false;
- int globalIdx = this->dofMapper().map(*eIt);
- const GlobalPosition &globalPos = eIt->geometry().center();
+ int globalIdx = this->dofMapper().map(*elemIt);
+ const GlobalPosition &globalPos = elemIt->geometry().center();
// initialize phase presence
staticDat_[globalIdx].phasePresence
@@ -215,10 +216,10 @@ public:
{
storage = 0;
- ElementIterator eIt = this->gridView_().template begin<0>();
- const ElementIterator eEndIt = this->gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ const ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ this->localResidual().evalPhaseStorage(*elemIt, phaseIdx);
for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
storage += this->localResidual().storageTerm()[i];
@@ -320,7 +321,7 @@ public:
VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
- if(velocityOutput_) // check if velocity output is demanded
+ if (velocityOutput_) // check if velocity output is demanded
{
// initialize velocity fields
for (unsigned int i = 0; i < numDofs; ++i)
@@ -337,25 +338,21 @@ public:
FVElementGeometry fvGeometry;
VolumeVariables volVars;
- ElementIterator eIt = this->gridView_().template begin<0>();
- ElementIterator eEndIt = this->gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt)
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt)
{
- int idx = this->elementMapper().map(*eIt);
+ int idx = this->elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
- fvGeometry.update(this->gridView_(), *eIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (numDofs == numElements) // element data
- globalIdx = idx;
- else
- globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
scvIdx,
false);
@@ -380,13 +377,13 @@ public:
(*temperature)[globalIdx] = volVars.temperature();
(*phasePresence)[globalIdx]
= staticDat_[globalIdx].phasePresence;
- if(velocityOutput_)
+ if (velocityOutput_)
{
(*cellNum)[globalIdx] += 1;
}
}
- if(velocityOutput_)
+ if (velocityOutput_)
{
// calculate vertex velocities
GlobalPosition tmpVelocity[numPhases];
@@ -405,7 +402,7 @@ public:
ElementVolumeVariables elemVolVars;
elemVolVars.update(this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
false /* oldSol? */);
@@ -413,7 +410,7 @@ public:
{
FluxVariables fluxVars(this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
faceIdx,
elemVolVars);
@@ -425,7 +422,7 @@ public:
// Transformation of the global normal vector to normal vector in the reference element
const typename Element::Geometry::JacobianTransposed jacobianT1 =
- eIt->geometry().jacobianTransposed(localPosIP);
+ elemIt->geometry().jacobianTransposed(localPosIP);
const GlobalPosition globalNormal = fluxVars.face().normal;
GlobalPosition localNormal(0);
@@ -444,11 +441,11 @@ public:
tmpVelocity[phaseIdx] = localNormal;
tmpVelocity[phaseIdx] *= flux[phaseIdx];
- if(phaseIdx == wPhaseIdx){
+ if (phaseIdx == wPhaseIdx){
scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
- else if(phaseIdx == nPhaseIdx){
+ else if (phaseIdx == nPhaseIdx){
scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
@@ -457,33 +454,33 @@ public:
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
const Dune::FieldVector<Scalar, dim>& localPos =
- ReferenceElements::general(eIt->geometry().type()).position(0, 0);
+ ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed jacobianT2 =
- eIt->geometry().jacobianTransposed(localPos);
+ elemIt->geometry().jacobianTransposed(localPos);
// transform vertex velocities from local to global coordinates
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocityW[scvIdx], scvVelocity);
- scvVelocity /= eIt->geometry().integrationElement(localPos);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the wetting phase subcontrolvolume velocities for each vertex
(*velocityW)[globalIdx] += scvVelocity;
jacobianT2.mtv(scvVelocityN[scvIdx], scvVelocity);
- scvVelocity /= eIt->geometry().integrationElement(localPos);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the nonwetting phase subcontrolvolume velocities for each vertex
(*velocityN)[globalIdx] += scvVelocity;
}
} // velocity output
} // loop over elements
- if(velocityOutput_)
+ if (velocityOutput_)
{
// divide the vertex velocities by the number of adjacent scvs i.e. cells
for(unsigned int globalIdx = 0; globalIdx < numDofs; ++globalIdx){
@@ -492,59 +489,32 @@ public:
}
}
- if (isBox) // vertex data
+ writer.attachDofData(*sN, "Sn", isBox);
+ writer.attachDofData(*sW, "Sw", isBox);
+ writer.attachDofData(*pN, "pN", isBox);
+ writer.attachDofData(*pW, "pW", isBox);
+ writer.attachDofData(*pC, "pC", isBox);
+ writer.attachDofData(*rhoW, "rhoW", isBox);
+ writer.attachDofData(*rhoN, "rhoN", isBox);
+ writer.attachDofData(*mobW, "mobW", isBox);
+ writer.attachDofData(*mobN, "mobN", isBox);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- 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");
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- std::ostringstream oss;
- oss << "X_" << FluidSystem::phaseName(phaseIdx) << "^" << FluidSystem::componentName(compIdx);
- writer.attachVertexData(*massFrac[phaseIdx][compIdx], oss.str());
- }
- }
- writer.attachVertexData(*poro, "porosity");
- writer.attachVertexData(*temperature, "temperature");
- writer.attachVertexData(*phasePresence, "phase presence");
-
- if(velocityOutput_) // check if velocity output is demanded
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
- writer.attachVertexData(*velocityW, "velocityW", dim);
- writer.attachVertexData(*velocityN, "velocityN", dim);
+ std::ostringstream oss;
+ oss << "X_" << FluidSystem::phaseName(phaseIdx) << "^" << FluidSystem::componentName(compIdx);
+ writer.attachDofData(*massFrac[phaseIdx][compIdx], oss.str(), isBox);
}
}
- else // cell data
+ writer.attachDofData(*poro, "porosity", isBox);
+ writer.attachDofData(*temperature, "temperature", isBox);
+ writer.attachDofData(*phasePresence, "phase presence", isBox);
+
+ if (velocityOutput_) // check if velocity output is demanded
{
- 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");
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- std::ostringstream oss;
- oss << "X_" << FluidSystem::phaseName(phaseIdx) << "^" << FluidSystem::componentName(compIdx);
- writer.attachCellData(*massFrac[phaseIdx][compIdx], oss.str());
- }
- }
- writer.attachCellData(*poro, "porosity");
- writer.attachCellData(*temperature, "temperature");
- writer.attachCellData(*phasePresence, "phase presence");
+ writer.attachDofData(*velocityW, "velocityW", isBox, dim);
+ writer.attachDofData(*velocityN, "velocityN", isBox, dim);
}
writer.attachCellData(*rank, "process rank");
@@ -609,19 +579,14 @@ public:
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
- ElementIterator eIt = this->gridView_().template begin<0> ();
- const ElementIterator &eEndIt = this->gridView_().template end<0> ();
- for (; eIt != eEndIt; ++eIt)
+ ElementIterator elemIt = this->gridView_().template begin<0> ();
+ const ElementIterator &elemEndIt = this->gridView_().template end<0> ();
+ for (; elemIt != elemEndIt; ++elemIt)
{
- fvGeometry.update(this->gridView_(), *eIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
-
- if (isBox)
- globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*eIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
if (staticDat_[globalIdx].visited)
continue;
@@ -629,7 +594,7 @@ public:
staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
scvIdx,
false);
diff --git a/dumux/implicit/2p2c/2p2cvolumevariables.hh b/dumux/implicit/2p2c/2p2cvolumevariables.hh
index 9591cae3f7..cc1d8bfd66 100644
--- a/dumux/implicit/2p2c/2p2cvolumevariables.hh
+++ b/dumux/implicit/2p2c/2p2cvolumevariables.hh
@@ -102,6 +102,7 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
//! The type of the object returned by the fluidState() method
@@ -184,14 +185,7 @@ public:
fvGeometry, scvIdx);
fluidState.setTemperature(t);
- unsigned numDofs = problem.model().numDofs();
- unsigned numVertices = problem.gridView().size(dim);
-
- int globalIdx;
- if (isBox) // vertex data
- globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
- else
- globalIdx = problem.model().dofMapper().map(element);
+ int globalIdx = problem.model().dofMapper().map(element, scvIdx, dofCodim);
int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
diff --git a/dumux/implicit/3p3c/3p3cmodel.hh b/dumux/implicit/3p3c/3p3cmodel.hh
index 9d8d48cd64..b973451b1c 100644
--- a/dumux/implicit/3p3c/3p3cmodel.hh
+++ b/dumux/implicit/3p3c/3p3cmodel.hh
@@ -140,6 +140,7 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -206,10 +207,10 @@ public:
{
storage = 0;
- ElementIterator eIt = this->gridView_().template begin<0>();
- const ElementIterator eEndIt = this->gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ const ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ this->localResidual().evalPhaseStorage(*elemIt, phaseIdx);
for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
storage += this->localResidual().storageTerm()[i];
@@ -325,11 +326,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox) // vertex data
- globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = idx;
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
@@ -362,64 +359,32 @@ public:
}
- if (isBox) // vertex data
- {
- writer.attachVertexData(*saturation[wPhaseIdx], "Sw");
- writer.attachVertexData(*saturation[nPhaseIdx], "Sn");
- writer.attachVertexData(*saturation[gPhaseIdx], "Sg");
- writer.attachVertexData(*pressure[wPhaseIdx], "pw");
- writer.attachVertexData(*pressure[nPhaseIdx], "pn");
- writer.attachVertexData(*pressure[gPhaseIdx], "pg");
- writer.attachVertexData(*density[wPhaseIdx], "rhow");
- writer.attachVertexData(*density[nPhaseIdx], "rhon");
- writer.attachVertexData(*density[gPhaseIdx], "rhog");
-
- for (int i = 0; i < numPhases; ++i)
- {
- for (int j = 0; j < numComponents; ++j)
- {
- std::ostringstream oss;
- oss << "x^"
- << FluidSystem::phaseName(i)
- << "_"
- << FluidSystem::componentName(j);
- writer.attachVertexData(*moleFraction[i][j], oss.str().c_str());
- }
- }
- writer.attachVertexData(*poro, "porosity");
- writer.attachVertexData(*perm, "permeability");
- writer.attachVertexData(*temperature, "temperature");
- writer.attachVertexData(*phasePresence, "phase presence");
- }
- else // cell data
+ writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
+ writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
+ writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);
+ writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
+ writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
+ writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
+ writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
+ writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
+ writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
+
+ for (int i = 0; i < numPhases; ++i)
{
- writer.attachCellData(*saturation[wPhaseIdx], "Sw");
- writer.attachCellData(*saturation[nPhaseIdx], "Sn");
- writer.attachCellData(*saturation[gPhaseIdx], "Sg");
- writer.attachCellData(*pressure[wPhaseIdx], "pw");
- writer.attachCellData(*pressure[nPhaseIdx], "pn");
- writer.attachCellData(*pressure[gPhaseIdx], "pg");
- writer.attachCellData(*density[wPhaseIdx], "rhow");
- writer.attachCellData(*density[nPhaseIdx], "rhon");
- writer.attachCellData(*density[gPhaseIdx], "rhog");
-
- for (int i = 0; i < numPhases; ++i)
- {
- for (int j = 0; j < numComponents; ++j)
- {
- std::ostringstream oss;
- oss << "x^"
- << FluidSystem::phaseName(i)
- << "_"
- << FluidSystem::componentName(j);
- writer.attachCellData(*moleFraction[i][j], oss.str().c_str());
- }
+ for (int j = 0; j < numComponents; ++j)
+ {
+ std::ostringstream oss;
+ oss << "x^"
+ << FluidSystem::phaseName(i)
+ << "_"
+ << FluidSystem::componentName(j);
+ writer.attachDofData(*moleFraction[i][j], oss.str().c_str(), isBox);
}
- writer.attachCellData(*poro, "porosity");
- writer.attachCellData(*perm, "permeability");
- writer.attachCellData(*temperature, "temperature");
- writer.attachCellData(*phasePresence, "phase presence");
}
+ writer.attachDofData(*poro, "porosity", isBox);
+ writer.attachDofData(*perm, "permeability", isBox);
+ writer.attachDofData(*temperature, "temperature", isBox);
+ writer.attachDofData(*phasePresence, "phase presence", isBox);
writer.attachCellData(*rank, "process rank");
}
@@ -483,19 +448,14 @@ public:
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
- ElementIterator eIt = this->gridView_().template begin<0> ();
- const ElementIterator &eEndIt = this->gridView_().template end<0> ();
- for (; eIt != eEndIt; ++eIt)
+ ElementIterator elemIt = this->gridView_().template begin<0> ();
+ const ElementIterator &elemEndIt = this->gridView_().template end<0> ();
+ for (; elemIt != elemEndIt; ++elemIt)
{
- fvGeometry.update(this->gridView_(), *eIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
-
- if (isBox)
- globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*eIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
if (staticDat_[globalIdx].visited)
continue;
@@ -503,7 +463,7 @@ public:
staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
scvIdx,
false);
diff --git a/dumux/implicit/3p3c/3p3cvolumevariables.hh b/dumux/implicit/3p3c/3p3cvolumevariables.hh
index edbbf9ebc8..52f0a55448 100644
--- a/dumux/implicit/3p3c/3p3cvolumevariables.hh
+++ b/dumux/implicit/3p3c/3p3cvolumevariables.hh
@@ -101,6 +101,7 @@ class ThreePThreeCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
static const Scalar R; // universial gas constant
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
//! The type of the object returned by the fluidState() method
@@ -130,11 +131,7 @@ public:
const MaterialLawParams &materialParams =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
- int globalIdx;
- if (isBox) // vertex data
- globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
- else
- globalIdx = problem.model().dofMapper().map(element);
+ int globalIdx = problem.model().dofMapper().map(element, scvIdx, dofCodim);
int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
diff --git a/dumux/implicit/co2/co2model.hh b/dumux/implicit/co2/co2model.hh
index c5f098d011..31b2fa649f 100644
--- a/dumux/implicit/co2/co2model.hh
+++ b/dumux/implicit/co2/co2model.hh
@@ -100,8 +100,10 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
- public:
+public:
/*!
@@ -118,24 +120,16 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
for (unsigned i = 0; i < ParentType::staticDat_.size(); ++i)
ParentType::staticDat_[i].visited = false;
- unsigned numDofs = this->numDofs();
- unsigned numVertices = this->problem_().gridView().size(dim);
-
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
- ElementIterator eIt = this->gridView_().template begin<0> ();
- const ElementIterator &eEndIt = this->gridView_().template end<0> ();
- for (; eIt != eEndIt; ++eIt)
+ ElementIterator elemIt = this->gridView_().template begin<0> ();
+ const ElementIterator &elemEndIt = this->gridView_().template end<0> ();
+ for (; elemIt != elemEndIt; ++elemIt)
{
- fvGeometry.update(this->gridView_(), *eIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
-
- if (numDofs != numVertices)
- globalIdx = this->elementMapper().map(*eIt);
- else
- globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
if (ParentType::staticDat_[globalIdx].visited)
continue;
@@ -143,17 +137,17 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
ParentType::staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
- *eIt,
+ *elemIt,
fvGeometry,
scvIdx,
false);
- const GlobalPosition &globalPos = eIt->geometry().corner(scvIdx);
+ const GlobalPosition &globalPos = elemIt->geometry().corner(scvIdx);
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
globalPos))
{
- this->jacobianAssembler().markVertexRed(globalIdx);
+ this->jacobianAssembler().markDofRed(globalIdx);
wasSwitched = true;
}
}
diff --git a/dumux/implicit/co2/co2volumevariables.hh b/dumux/implicit/co2/co2volumevariables.hh
index 3d346443ec..e494aa6510 100644
--- a/dumux/implicit/co2/co2volumevariables.hh
+++ b/dumux/implicit/co2/co2volumevariables.hh
@@ -87,6 +87,8 @@ class CO2VolumeVariables: public TwoPTwoCVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dim = GridView::dimension};
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
static const Scalar R; // universial nonwetting constant
@@ -112,7 +114,7 @@ public:
const int scvIdx,
const bool isOldSol)
{
- // Update BoxVolVars but not 2p2cvolvars
+ // Update BoxVolVars but not 2p2cvolvars
// ToDo: Is BaseClassType the right name?
BaseClassType::update(priVars,
problem,
@@ -121,14 +123,7 @@ public:
scvIdx,
isOldSol);
- unsigned numDofs = problem.model().numDofs();
- unsigned numVertices = problem.gridView().size(dim);
-
- int globalIdx;
- if (numDofs != numVertices) // element data
- globalIdx = problem.model().dofMapper().map(element);
- else
- globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ int globalIdx = problem.model().dofMapper().map(element, scvIdx, dofCodim);
int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
diff --git a/dumux/implicit/common/implicitmodel.hh b/dumux/implicit/common/implicitmodel.hh
index ceafafa3d1..1b45cb14e9 100644
--- a/dumux/implicit/common/implicitmodel.hh
+++ b/dumux/implicit/common/implicitmodel.hh
@@ -78,6 +78,7 @@ class ImplicitModel
typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
// copying a model is not a good idea
ImplicitModel(const ImplicitModel &);
@@ -884,18 +885,8 @@ protected:
// loop over all element vertices, i.e. sub control volumes
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; scvIdx++)
{
- int globalIdx;
-
- if (isBox)
- {
- // map the local vertex index to the global one
- globalIdx = vertexMapper().map(*eIt, scvIdx, dim);
- }
- else
- {
- // get the global index of the element
- globalIdx = elementMapper().map(*eIt);
- }
+ // get the global index of the degree of freedom
+ int globalIdx = dofMapper().map(*eIt, scvIdx, dofCodim);
// let the problem do the dirty work of nailing down
// the initial solution.
diff --git a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
index 79be576f16..ab71e76017 100644
--- a/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
+++ b/dumux/implicit/mpnc/energy/mpncvtkwriterenergy.hh
@@ -62,6 +62,7 @@ class MPNCVtkWriterEnergy : public MPNCVtkWriterModule<TypeTag>
typedef typename ParentType::ScalarVector ScalarVector;
typedef typename ParentType::PhaseVector PhaseVector;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
MPNCVtkWriterEnergy(const Problem &problem)
@@ -84,15 +85,14 @@ public:
* \brief Modify the internal buffers according to the volume
* variables seen on an element
*/
- void processElement(const Element &elem,
+ void processElement(const Element &element,
const FVElementGeometry &fvGeometry,
const ElementVolumeVariables &elemVolVars,
const ElementBoundaryTypes &elemBcTypes)
{
- int numLocalVertices = elem.geometry().corners();
- for (int localVertexIdx = 0; localVertexIdx < numLocalVertices; ++localVertexIdx) {
- const unsigned int globalIdx = this->problem_.vertexMapper().map(elem, localVertexIdx, dim);
- const VolumeVariables &volVars = elemVolVars[localVertexIdx];
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx) {
+ const unsigned int globalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
if (temperatureOutput_)
temperature_[globalIdx] = volVars.fluidState().temperature(/*phaseIdx=*/0);
@@ -144,6 +144,7 @@ class MPNCVtkWriterEnergy<TypeTag, /* enableEnergy = */ true, /* enableKineticEn
typedef typename ParentType::ScalarVector ScalarVector;
typedef typename ParentType::PhaseVector PhaseVector;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
MPNCVtkWriterEnergy(const Problem &problem)
@@ -170,15 +171,14 @@ public:
* \brief Modify the internal buffers according to the volume
* variables seen on an element
*/
- void processElement(const Element &elem,
+ void processElement(const Element &element,
const FVElementGeometry &fvGeometry,
const ElementVolumeVariables &elemVolVars,
const ElementBoundaryTypes &elemBcTypes)
{
- const unsigned int numLocalVertices = elem.geometry().corners();
- for (int localVertexIdx = 0; localVertexIdx < numLocalVertices; ++localVertexIdx) {
- int gobalIdx = this->problem_.vertexMapper().map(elem, localVertexIdx, dim);
- const VolumeVariables &volVars = elemVolVars[localVertexIdx];
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx) {
+ int gobalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
if (temperatureOutput_) temperature_[gobalIdx] = volVars.fluidState().temperature(/*phaseIdx=*/0);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
diff --git a/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh b/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh
index 3968f741df..e20077b6d3 100644
--- a/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh
+++ b/dumux/implicit/mpnc/mass/mpncvtkwritermass.hh
@@ -57,6 +57,7 @@ class MPNCVtkWriterMass : public MPNCVtkWriterModule<TypeTag>
enum { dim = GridView::dimension };
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
typedef typename ParentType::ComponentVector ComponentVector;
bool fugacityOutput_;
@@ -82,15 +83,14 @@ public:
* \brief Modify the internal buffers according to the volume
* variables seen on an element
*/
- void processElement(const Element &elem,
- const FVElementGeometry &fvElemGeom,
+ void processElement(const Element &element,
+ const FVElementGeometry &fvGeometry,
const ElementVolumeVariables &elemVolVars,
const ElementBoundaryTypes &elemBcTypes)
{
- int numLocalVertices = elem.geometry().corners();
- for (int localVertexIdx = 0; localVertexIdx < numLocalVertices; ++localVertexIdx) {
- int globalIdx = this->problem_.vertexMapper().map(elem, localVertexIdx, dim);
- const VolumeVariables &volVars = elemVolVars[localVertexIdx];
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx) {
+ const unsigned int globalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
if (fugacityOutput_) {
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
diff --git a/dumux/implicit/mpnc/mpncvtkwritercommon.hh b/dumux/implicit/mpnc/mpncvtkwritercommon.hh
index c6e6f74e10..703c9c8b75 100644
--- a/dumux/implicit/mpnc/mpncvtkwritercommon.hh
+++ b/dumux/implicit/mpnc/mpncvtkwritercommon.hh
@@ -63,6 +63,7 @@ class MPNCVtkWriterCommon : public MPNCVtkWriterModule<TypeTag>
typedef Dune::BlockVector<DimVector> DimField;
typedef Dune::array<DimField, numPhases> PhaseDimField;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
MPNCVtkWriterCommon(const Problem &problem)
@@ -127,11 +128,7 @@ public:
{
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox) // vertex data
- globalIdx = this->problem_.vertexMapper().map(element, scvIdx, dim);
- else
- globalIdx = this->problem_.elementMapper().map(element);
+ int globalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
const VolumeVariables &volVars = elemVolVars[scvIdx];
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index 8ed12767a9..7551403fc7 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -120,6 +120,7 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
enum { dim = GridView::dimension };
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -168,11 +169,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox)
- globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = this->elementMapper().map(*elemIt);
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
volVars.update(sol[globalIdx],
this->problem_(),
@@ -195,34 +192,17 @@ public:
}
}
- if (isBox) // 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");
- }
- else // cell 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");
- }
+ writer.attachDofData(*Sn, "Sn", isBox);
+ writer.attachDofData(*Sw, "Sw", isBox);
+ writer.attachDofData(*pN, "pn", isBox);
+ writer.attachDofData(*pW, "pw", isBox);
+ writer.attachDofData(*pC, "pc", isBox);
+ writer.attachDofData(*rhoW, "rhoW", isBox);
+ writer.attachDofData(*rhoN, "rhoN", isBox);
+ writer.attachDofData(*mobW, "mobW", isBox);
+ writer.attachDofData(*mobN, "mobN", isBox);
+ writer.attachDofData(*poro, "porosity", isBox);
+ writer.attachDofData(*Te, "temperature", isBox);
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/richards/richardsnewtoncontroller.hh b/dumux/implicit/richards/richardsnewtoncontroller.hh
index 2f0b035d6d..5513dcf5af 100644
--- a/dumux/implicit/richards/richardsnewtoncontroller.hh
+++ b/dumux/implicit/richards/richardsnewtoncontroller.hh
@@ -55,6 +55,7 @@ class RichardsNewtonController : public NewtonController<TypeTag>
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum { dim = GridView::dimension };
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
/*!
@@ -96,11 +97,7 @@ public:
fvGeometry.update(gridView, *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx;
- if (isBox)
- globalIdx = this->problem_().vertexMapper().map(*elemIt, scvIdx, dim);
- else
- globalIdx = this->problem_().elementMapper().map(*elemIt);
+ int globalIdx = this->model_().dofMapper().map(*elemIt, scvIdx, dofCodim);
// calculate the old wetting phase saturation
const SpatialParams &spatialParams = this->problem_().spatialParams();
diff --git a/dumux/io/vtkmultiwriter.hh b/dumux/io/vtkmultiwriter.hh
index c4c65b0c80..00c0cb0149 100644
--- a/dumux/io/vtkmultiwriter.hh
+++ b/dumux/io/vtkmultiwriter.hh
@@ -219,6 +219,22 @@ public:
curWriter_->addCellData(fnPtr);
}
+ /*!
+ * \brief Add data associated with degrees of freedom to the output.
+ *
+ * This is a wrapper for the functions attachVertexData and attachCelldata.
+ * Depending on the value of \param isVertexData, the appropriate function
+ * is selected.
+ */
+ template <class DataBuffer>
+ void attachDofData(DataBuffer &buf, std::string name, bool isVertexData, int nComps = 1)
+ {
+ if (isVertexData)
+ attachVertexData(buf, name, nComps);
+ else
+ attachCellData(buf, name, nComps);
+ }
+
/*!
* \brief Finalizes the current writer.
*
--
GitLab