Commit bb0e3a30 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

implicit branch: add attachDofData method to the VTKMultiWriter to uniformly...

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
parent 764429db
......@@ -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");
}
};
......
......@@ -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");
}
......
......@@ -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");
}
......
......@@ -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: