diff --git a/dumux/implicit/2p/2pmodel.hh b/dumux/implicit/2p/2pmodel.hh
index 5a46820d2c3375556e9af53b3978061145af7214..40922141833b1ec0cc674b392341c20f404fe972 100644
--- a/dumux/implicit/2p/2pmodel.hh
+++ b/dumux/implicit/2p/2pmodel.hh
@@ -98,6 +98,8 @@ class TwoPModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief Append all quantities of interest which can be derived
@@ -115,14 +117,13 @@ public:
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
- // get the number of degrees of freedom and the number of vertices
+ // get the number of degrees of freedom
unsigned numDofs = this->numDofs();
- unsigned numVertices = this->problem_().gridView().size(dim);
-
- // velocity output currently only works for vertex data
- if (numDofs != numVertices)
- velocityOutput = false;
-
+
+ // velocity output currently only works for the box discretization
+ if (!isBox)
+ velocityOutput = false;
+
// create the required scalar fields
ScalarField *pW = writer.allocateManagedBuffer(numDofs);
ScalarField *pN = writer.allocateManagedBuffer(numDofs);
@@ -142,7 +143,7 @@ public:
if(velocityOutput) // check if velocity output is demanded
{
// initialize velocity fields
- for (unsigned int i = 0; i < numVertices; ++i)
+ for (unsigned int i = 0; i < numDofs; ++i)
{
(*velocityN)[i] = Scalar(0);
(*velocityW)[i] = Scalar(0);
@@ -170,169 +171,135 @@ public:
fvGeometry.update(this->gridView_(), *elemIt);
- if (numDofs == numElements) // element data
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- volVars.update(sol[idx],
+ int globalIdx;
+ if (isBox)
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*elemIt);
+
+ volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
fvGeometry,
- /*scvIdx=*/0,
+ scvIdx,
false);
- (*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();
+ (*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;
+ }
}
- else // vertex data
+
+ if(velocityOutput)
{
- int numVerts = elemIt->template count<dim> ();
- for (int scvIdx = 0; scvIdx < numVerts; ++scvIdx)
+ // calculate vertex velocities
+ GlobalPosition tmpVelocity[numPhases];
+
+ for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
- volVars.update(sol[globalIdx],
- this->problem_(),
- *elemIt,
- fvGeometry,
- scvIdx,
- 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;
- }
+ tmpVelocity[phaseIdx] = Scalar(0.0);
}
- if(velocityOutput)
- {
- // calculate vertex velocities
- GlobalPosition tmpVelocity[numPhases];
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
+ SCVVelocities scvVelocityW(8), scvVelocityN(8);
- for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- tmpVelocity[phaseIdx] = Scalar(0.0);
- }
+ scvVelocityW = 0;
+ scvVelocityN = 0;
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
- SCVVelocities scvVelocityW(8), scvVelocityN(8);
+ ElementVolumeVariables elemVolVars;
- scvVelocityW = 0;
- scvVelocityN = 0;
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ false /* oldSol? */);
- ElementVolumeVariables elemVolVars;
+ for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
+ {
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvGeometry,
- false /* oldSol? */);
+ FluxVariables fluxVars(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ faceIdx,
+ elemVolVars);
- for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
-
- FluxVariables fluxVars(this->problem_(),
- *elemIt,
- fvGeometry,
- faceIdx,
- elemVolVars);
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- // local position of integration point
- const Dune::FieldVector<Scalar, dim>& localPosIP = fvGeometry.subContVolFace[faceIdx].ipLocal;
-
- // Transformation of the global normal vector to normal vector in the reference element
- const typename Element::Geometry::JacobianTransposed& jacobianT1 =
- elemIt->geometry().jacobianTransposed(localPosIP);
-
- const GlobalPosition globalNormal = fluxVars.face().normal;
- GlobalPosition localNormal;
- 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.
- PhasesVector q;
- q[phaseIdx] = fluxVars.volumeFlux(phaseIdx) / localArea;
-
- // transform the normal Darcy velocity into a vector
- tmpVelocity[phaseIdx] = localNormal;
- tmpVelocity[phaseIdx] *= q[phaseIdx];
-
- if(phaseIdx == wPhaseIdx){
- scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
- scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
- }
- else if(phaseIdx == nPhaseIdx){
- scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
- scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
- }
- }
+ // local position of integration point
+ const Dune::FieldVector<Scalar, dim>& localPosIP = fvGeometry.subContVolFace[faceIdx].ipLocal;
+
+ // Transformation of the global normal vector to normal vector in the reference element
+ const typename Element::Geometry::JacobianTransposed& jacobianT1 =
+ elemIt->geometry().jacobianTransposed(localPosIP);
+
+ const GlobalPosition globalNormal = fluxVars.face().normal;
+ GlobalPosition localNormal;
+ 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.
+ PhasesVector q;
+ q[phaseIdx] = fluxVars.volumeFlux(phaseIdx) / localArea;
+
+ // transform the normal Darcy velocity into a vector
+ tmpVelocity[phaseIdx] = localNormal;
+ tmpVelocity[phaseIdx] *= q[phaseIdx];
+
+ if(phaseIdx == wPhaseIdx){
+ scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
+ }
+ else if(phaseIdx == nPhaseIdx){
+ scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
+ }
}
- typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
- const Dune::FieldVector<Scalar, dim> &localPos
- = ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
+ }
+ 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 typename Element::Geometry::JacobianTransposed& jacobianT2
- = elemIt->geometry().jacobianTransposed(localPos);
+ // get the transposed Jacobian of the element mapping
+ const typename Element::Geometry::JacobianTransposed& 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 vertex velocities from local to global coordinates
+ for (int i = 0; i < fvGeometry.numSCV; ++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 (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
+ if (isBox) // vertex data
{
writer.attachVertexData(*Sn, "Sn");
writer.attachVertexData(*Sw, "Sw");
@@ -348,7 +315,7 @@ public:
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 < numVertices; ++globalIdx){
+ for(unsigned int globalIdx = 0; globalIdx < velocityW->size(); ++globalIdx){
(*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
(*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
}
@@ -356,6 +323,20 @@ public:
writer.attachVertexData(*velocityN, "velocityN", dim);
}
}
+ 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.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index b3b395a54dc607c5bf8ebc0f66ec10fefe70ce61..82a91f9e163adf24780ba9419cebdc846fc84550 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -134,6 +134,8 @@ class TwoPTwoCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief Initialize the static data with the initial solution.
@@ -145,9 +147,8 @@ public:
ParentType::init(problem);
unsigned numDofs = this->numDofs();
- unsigned numVertices = this->problem_().gridView().size(dim);
- staticVertexDat_.resize(numDofs);
+ staticDat_.resize(numDofs);
setSwitched_(false);
@@ -161,7 +162,7 @@ public:
velocityOutput_ = false;
}
- if (numDofs != numVertices) // i.e. cell-centered discretization
+ if (!isBox) // i.e. cell-centered discretization
{
velocityOutput_ = false;
@@ -169,20 +170,20 @@ public:
const GlobalPosition &globalPos = eIt->geometry().center();
// initialize phase presence
- staticVertexDat_[globalIdx].phasePresence
+ staticDat_[globalIdx].phasePresence
= this->problem_().initialPhasePresence(*(this->gridView_().template begin<dim>()),
globalIdx, globalPos);
- staticVertexDat_[globalIdx].wasSwitched = false;
+ staticDat_[globalIdx].wasSwitched = false;
- staticVertexDat_[globalIdx].oldPhasePresence
- = staticVertexDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
}
}
if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
- if (numDofs == numVertices) // i.e. vertex-centered discretization
+ if (isBox) // i.e. vertex-centered discretization
{
VertexIterator vIt = this->gridView_().template begin<dim> ();
const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
@@ -192,13 +193,13 @@ public:
const GlobalPosition &globalPos = vIt->geometry().corner(0);
// initialize phase presence
- staticVertexDat_[globalIdx].phasePresence
+ staticDat_[globalIdx].phasePresence
= this->problem_().initialPhasePresence(*vIt, globalIdx,
globalPos);
- staticVertexDat_[globalIdx].wasSwitched = false;
+ staticDat_[globalIdx].wasSwitched = false;
- staticVertexDat_[globalIdx].oldPhasePresence
- = staticVertexDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
}
}
}
@@ -272,8 +273,8 @@ public:
*/
int phasePresence(int globalIdx, bool oldSol) const
{
- return oldSol ? staticVertexDat_[globalIdx].oldPhasePresence
- : staticVertexDat_[globalIdx].phasePresence;
+ return oldSol ? staticDat_[globalIdx].oldPhasePresence
+ : staticDat_[globalIdx].phasePresence;
}
/*!
@@ -291,14 +292,14 @@ 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
unsigned numDofs = this->numDofs();
- unsigned numVertices = this->problem_().gridView().size(dim);
- // velocity output currently only works for vertex data
- if (numDofs != numVertices)
+ // velocity output currently only works for the box discretization
+ if (!isBox)
velocityOutput_ = false;
+ // create the required scalar fields
ScalarField *sN = writer.allocateManagedBuffer(numDofs);
ScalarField *sW = writer.allocateManagedBuffer(numDofs);
ScalarField *pN = writer.allocateManagedBuffer(numDofs);
@@ -378,7 +379,7 @@ public:
(*poro)[globalIdx] = volVars.porosity();
(*temperature)[globalIdx] = volVars.temperature();
(*phasePresence)[globalIdx]
- = staticVertexDat_[globalIdx].phasePresence;
+ = staticDat_[globalIdx].phasePresence;
if(velocityOutput_)
{
(*cellNum)[globalIdx] += 1;
@@ -491,31 +492,7 @@ public:
}
}
- 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");
- 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");
- }
- else
+ if (isBox) // vertex data
{
writer.attachVertexData(*sN, "Sn");
writer.attachVertexData(*sW, "Sw");
@@ -545,6 +522,30 @@ public:
writer.attachVertexData(*velocityN, "velocityN", dim);
}
}
+ 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");
+ 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.attachCellData(*rank, "process rank");
}
@@ -565,7 +566,7 @@ public:
if (!outStream.good())
DUNE_THROW(Dune::IOError, "Could not serialize entity " << globalIdx);
- outStream << staticVertexDat_[globalIdx].phasePresence << " ";
+ outStream << staticDat_[globalIdx].phasePresence << " ";
}
/*!
@@ -586,9 +587,9 @@ public:
DUNE_THROW(Dune::IOError,
"Could not deserialize entity " << globalIdx);
- inStream >> staticVertexDat_[globalIdx].phasePresence;
- staticVertexDat_[globalIdx].oldPhasePresence
- = staticVertexDat_[globalIdx].phasePresence;
+ inStream >> staticDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
}
@@ -603,11 +604,8 @@ public:
{
bool wasSwitched = false;
- for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
- staticVertexDat_[i].visited = false;
-
- unsigned numDofs = this->numDofs();
- unsigned numVertices = this->problem_().gridView().size(dim);
+ for (unsigned i = 0; i < staticDat_.size(); ++i)
+ staticDat_[i].visited = false;
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
@@ -620,28 +618,28 @@ public:
{
int globalIdx;
- if (numDofs != numVertices)
- globalIdx = this->elementMapper().map(*eIt);
- else
+ if (isBox)
globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*eIt);
- if (staticVertexDat_[globalIdx].visited)
+ if (staticDat_[globalIdx].visited)
continue;
- staticVertexDat_[globalIdx].visited = true;
+ staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
*eIt,
fvGeometry,
scvIdx,
false);
- const GlobalPosition &globalPos = eIt->geometry().corner(scvIdx);
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
globalPos))
{
- this->jacobianAssembler().markVertexRed(globalIdx);
+ this->jacobianAssembler().markDofRed(globalIdx);
wasSwitched = true;
}
}
@@ -677,11 +675,11 @@ protected:
*/
void resetPhasePresence_()
{
- for (unsigned int idx = 0; idx < staticVertexDat_.size(); ++idx)
+ for (unsigned int idx = 0; idx < staticDat_.size(); ++idx)
{
- staticVertexDat_[idx].phasePresence
- = staticVertexDat_[idx].oldPhasePresence;
- staticVertexDat_[idx].wasSwitched = false;
+ staticDat_[idx].phasePresence
+ = staticDat_[idx].oldPhasePresence;
+ staticDat_[idx].wasSwitched = false;
}
}
@@ -690,11 +688,11 @@ protected:
*/
void updateOldPhasePresence_()
{
- for (unsigned int idx = 0; idx < staticVertexDat_.size(); ++idx)
+ for (unsigned int idx = 0; idx < staticDat_.size(); ++idx)
{
- staticVertexDat_[idx].oldPhasePresence
- = staticVertexDat_[idx].phasePresence;
- staticVertexDat_[idx].wasSwitched = false;
+ staticDat_[idx].oldPhasePresence
+ = staticDat_[idx].phasePresence;
+ staticDat_[idx].wasSwitched = false;
}
}
@@ -715,7 +713,7 @@ protected:
{
// evaluate primary variable switch
bool wouldSwitch = false;
- int phasePresence = staticVertexDat_[globalIdx].phasePresence;
+ int phasePresence = staticDat_[globalIdx].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
@@ -728,7 +726,7 @@ protected:
Scalar xwMax = 1.0;
if (xww + xwn > xwMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xwMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -756,7 +754,7 @@ protected:
Scalar xgMax = 1.0;
if (xnw + xnn > xgMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xgMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -777,7 +775,7 @@ protected:
else if (phasePresence == bothPhases)
{
Scalar Smin = 0.0;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
@@ -806,14 +804,14 @@ protected:
}
}
- staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
- staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
+ staticDat_[globalIdx].phasePresence = newPhasePresence;
+ staticDat_[globalIdx].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
protected:
// parameters given in constructor
- std::vector<StaticVars> staticVertexDat_;
+ std::vector<StaticVars> staticDat_;
bool switchFlag_;
bool velocityOutput_;
};
diff --git a/dumux/implicit/co2/co2model.hh b/dumux/implicit/co2/co2model.hh
index 1d3f2edb82b92a54af3c230a5290f0c70c6ccc3c..c5f098d011eaaea27d9eaf9ad9ed5f4b255d5e8d 100644
--- a/dumux/implicit/co2/co2model.hh
+++ b/dumux/implicit/co2/co2model.hh
@@ -115,8 +115,8 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
{
bool wasSwitched = false;
- for (unsigned i = 0; i < ParentType::staticVertexDat_.size(); ++i)
- ParentType::staticVertexDat_[i].visited = false;
+ 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);
@@ -137,10 +137,10 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
else
globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
- if (ParentType::staticVertexDat_[globalIdx].visited)
+ if (ParentType::staticDat_[globalIdx].visited)
continue;
- ParentType::staticVertexDat_[globalIdx].visited = true;
+ ParentType::staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
*eIt,
@@ -181,7 +181,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
typename FluidSystem::ParameterCache paramCache;
// evaluate primary variable switch
bool wouldSwitch = false;
- int phasePresence = ParentType::staticVertexDat_[globalIdx].phasePresence;
+ int phasePresence = ParentType::staticDat_[globalIdx].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
@@ -194,7 +194,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
if(xnw > xnwMax)
wouldSwitch = true;
- if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
+ if (ParentType::staticDat_[globalIdx].wasSwitched)
xnwMax *= 1.02;
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
@@ -220,7 +220,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if(xwn > xwnMax)
wouldSwitch = true;
- if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
+ if (ParentType::staticDat_[globalIdx].wasSwitched)
xwnMax *= 1.02;
@@ -241,7 +241,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
else if (phasePresence == bothPhases)
{
Scalar Smin = 0.0;
- if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
+ if (ParentType::staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
@@ -270,8 +270,8 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
}
}
- ParentType::staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
- ParentType::staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
+ ParentType::staticDat_[globalIdx].phasePresence = newPhasePresence;
+ ParentType::staticDat_[globalIdx].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
diff --git a/dumux/implicit/common/implicitassembler.hh b/dumux/implicit/common/implicitassembler.hh
index 9c731fe616cfe78a99968b874fba0bf3bb646d24..218cae454c9bd50f6578ffc7237f70b5877e437d 100644
--- a/dumux/implicit/common/implicitassembler.hh
+++ b/dumux/implicit/common/implicitassembler.hh
@@ -298,6 +298,23 @@ public:
}
+ /*!
+ * \brief Force to reassemble a given degree of freedom
+ * next time the assemble() method is called.
+ *
+ * \param globalIdx The global index of the degree of freedom
+ */
+ void markDofRed(const int globalIdx)
+ {
+ if (!enablePartialReassemble_())
+ return;
+
+ if (isBox)
+ vertexColor_[globalIdx] = Red;
+ else
+ elementColor_[globalIdx] = Red;
+ }
+
/*!
* \brief Force to reassemble a given vertex next time the
* assemble() method is called.
@@ -305,6 +322,7 @@ public:
* \param globalVertIdx The global index of the vertex which ought
* to be red.
*/
+ DUNE_DEPRECATED_MSG("Use markDofRed instead.")
void markVertexRed(const int globalVertIdx)
{
if (!enablePartialReassemble_())