diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh index c3e229974bff92ab3b2b5cd8f681e3f02fa95f67..b51cbae2f5a65b373abe76253e2528deda1dd521 100644 --- a/dumux/boxmodels/common/boxelementboundarytypes.hh +++ b/dumux/boxmodels/common/boxelementboundarytypes.hh @@ -81,11 +81,9 @@ public: * \param problem The problem object which needs to be simulated * \param element The DUNE Codim<0> entity for which the boundary * types should be collected - * \param fvElemGeom The element's finite volume geometry */ void update(const Problem &problem, - const Element &element, - const FVElementGeometry &fvElemGeom) + const Element &element) { int numVerts = element.template count<dim>(); this->resize(numVerts); @@ -109,6 +107,20 @@ public: } } + /*! + * \brief Update the boundary types for all vertices of an element. + * + * \param problem The problem object which needs to be simulated + * \param element The DUNE Codim<0> entity for which the boundary + * types should be collected + * \param fvGeometry The element's finite volume geometry + */ + DUMUX_DEPRECATED_MSG("use update(problem, element) instead") + void update(const Problem &problem, + const Element &element, + const FVElementGeometry &fvGeometry) + { update(problem, element); } + /*! * \brief Returns whether the element has a vertex which contains * a Dirichlet value. diff --git a/dumux/boxmodels/common/boxelementvolumevariables.hh b/dumux/boxmodels/common/boxelementvolumevariables.hh index bed72fa53111d5240c71c0d5056ec781a21d562b..31d02d29fe5209b2ec380d1148af54141fe9c6f6 100644 --- a/dumux/boxmodels/common/boxelementvolumevariables.hh +++ b/dumux/boxmodels/common/boxelementvolumevariables.hh @@ -66,12 +66,12 @@ public: * * \param problem The problem which needs to be simulated. * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated - * \param fvElemGeom The finite volume geometry of the element + * \param fvGeometry The finite volume geometry of the element * \param oldSol Tells whether the model's previous or current solution should be used. */ void update(const Problem &problem, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, bool oldSol) { const SolutionVector &globalSol = @@ -81,16 +81,16 @@ public: const VertexMapper &vertexMapper = problem.vertexMapper(); // we assert that the i-th shape function is // associated to the i-th vert of the element. - int n = element.template count<dim>(); - this->resize(n); - for (int i = 0; i < n; i++) { + int numVertices = element.template count<dim>(); + this->resize(numVertices); + for (int scvIdx = 0; scvIdx < numVertices; scvIdx++) { const PrimaryVariables &solI - = globalSol[vertexMapper.map(element, i, dim)]; - (*this)[i].update(solI, + = globalSol[vertexMapper.map(element, scvIdx, dim)]; + (*this)[scvIdx].update(solI, problem, element, - fvElemGeom, - i, + fvGeometry, + scvIdx, oldSol); } } @@ -99,33 +99,33 @@ public: * \brief Construct the volume variables for all of vertices of an * element given a solution vector computed by PDELab. * - * \tparam ElemSolVectorType The container type which stores the + * \tparam ElementSolutionVector The container type which stores the * primary variables of the element * using _local_ indices * * \param problem The problem which needs to be simulated. * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated - * \param fvElemGeom The finite volume geometry of the element + * \param fvGeometry The finite volume geometry of the element * \param elementSolVector The local solution for the element using PDELab ordering */ - template<typename ElemSolVectorType> + template<typename ElementSolutionVector> void updatePDELab(const Problem &problem, const Element &element, - const FVElementGeometry &fvElemGeom, - const ElemSolVectorType& elementSolVector) + const FVElementGeometry &fvGeometry, + const ElementSolutionVector& elementSolVector) { - int n = element.template count<dim>(); - this->resize(n); - for (int vertexIdx = 0; vertexIdx < n; vertexIdx++) + int numVertices = element.template count<dim>(); + this->resize(numVertices); + for (int scvIdx = 0; scvIdx < numVertices; scvIdx++) { PrimaryVariables solI(0); - for (int eqnIdx=0; eqnIdx<numEq; eqnIdx++) - solI[eqnIdx] = elementSolVector[vertexIdx + eqnIdx*n]; - (*this)[vertexIdx].update(solI, + for (int eqnIdx = 0; eqnIdx < numEq; eqnIdx++) + solI[eqnIdx] = elementSolVector[scvIdx + eqnIdx*numVertices]; + (*this)[scvIdx].update(solI, problem, element, - fvElemGeom, - vertexIdx, + fvGeometry, + scvIdx, false); } diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh index c856138b22d2235da4a5b34eb1dcf315ffe959bd..7760bfd5b4f50f0276e05eb513149497dadd4fda 100644 --- a/dumux/boxmodels/common/boxlocaljacobian.hh +++ b/dumux/boxmodels/common/boxlocaljacobian.hh @@ -141,14 +141,14 @@ public: fvElemGeom_.update(gridView_(), element); reset_(); - bcTypes_.update(problem_(), elem_(), fvElemGeom_); + bcTypes_.update(problem_(), element_(), fvElemGeom_); // this is pretty much a HACK because the internal state of // the problem is not supposed to be changed during the // evaluation of the residual. (Reasons: It is a violation of // abstraction, makes everything more prone to errors and is // not thread save.) The real solution are context objects! - problem_().updateCouplingParams(elem_()); + problem_().updateCouplingParams(element_()); // set the hints for the volume variables model_().setHints(element, prevVolVars_, curVolVars_); @@ -156,12 +156,12 @@ public: // update the secondary variables for the element at the last // and the current time levels prevVolVars_.update(problem_(), - elem_(), + element_(), fvElemGeom_, true /* isOldSol? */); curVolVars_.update(problem_(), - elem_(), + element_(), fvElemGeom_, false /* isOldSol? */); @@ -169,7 +169,7 @@ public: model_().updateCurHints(element, curVolVars_); // calculate the local residual - localResidual().eval(elem_(), + localResidual().eval(element_(), fvElemGeom_, prevVolVars_, curVolVars_, @@ -177,7 +177,7 @@ public: residual_ = localResidual().residual(); storageTerm_ = localResidual().storageTerm(); - model_().updatePVWeights(elem_(), curVolVars_); + model_().updatePVWeights(element_(), curVolVars_); // calculate the local jacobian matrix int numVertices = fvElemGeom_.numVertices; @@ -303,6 +303,16 @@ protected: /*! * \brief Returns a reference to the element. */ + const Element &element_() const + { + Valgrind::CheckDefined(elemPtr_); + return *elemPtr_; + }; + + /*! + * \brief Returns a reference to the element. + */ + DUMUX_DEPRECATED_MSG("use element_() instead") const Element &elem_() const { Valgrind::CheckDefined(elemPtr_); @@ -332,7 +342,7 @@ protected: */ void reset_() { - int n = elem_().template count<dim>(); + int n = element_().template count<dim>(); for (int i = 0; i < n; ++ i) { storageJacobian_[i] = 0.0; for (int j = 0; j < n; ++ j) { @@ -375,9 +385,9 @@ protected: * is the value of a sub-control volume's primary variable at the * evaluation point and \f$\epsilon\f$ is a small value larger than 0. * - * \param dest The vector storing the partial derivatives of all + * \param partialDeriv The vector storing the partial derivatives of all * equations - * \param destStorage the mass matrix contributions + * \param storageDeriv the mass matrix contributions * \param scvIdx The sub-control volume index of the current * finite element for which the partial derivative * ought to be calculated @@ -386,12 +396,12 @@ protected: * for which the partial derivative ought to be * calculated */ - void evalPartialDerivative_(ElementSolutionVector &dest, - PrimaryVariables &destStorage, + void evalPartialDerivative_(ElementSolutionVector &partialDeriv, + PrimaryVariables &storageDeriv, int scvIdx, int pvIdx) { - int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim); + int globalIdx = vertexMapper_().map(element_(), scvIdx, dim); PrimaryVariables priVars(model_().curSol()[globalIdx]); VolumeVariables origVolVars(curVolVars_[scvIdx]); @@ -411,26 +421,26 @@ protected: // calculate the residual curVolVars_[scvIdx].update(priVars, problem_(), - elem_(), + element_(), fvElemGeom_, scvIdx, false); - localResidual().eval(elem_(), + localResidual().eval(element_(), fvElemGeom_, prevVolVars_, curVolVars_, bcTypes_); // store the residual and the storage term - dest = localResidual().residual(); - destStorage = localResidual().storageTerm()[scvIdx]; + partialDeriv = localResidual().residual(); + storageDeriv = localResidual().storageTerm()[scvIdx]; } else { // we are using backward differences, i.e. we don't need // to calculate f(x + \epsilon) and we can recycle the // (already calculated) residual f(x) - dest = residual_; - destStorage = storageTerm_[scvIdx]; + partialDeriv = residual_; + storageDeriv = storageTerm_[scvIdx]; } @@ -445,37 +455,37 @@ protected: // calculate residual again curVolVars_[scvIdx].update(priVars, problem_(), - elem_(), + element_(), fvElemGeom_, scvIdx, false); - localResidual().eval(elem_(), + localResidual().eval(element_(), fvElemGeom_, prevVolVars_, curVolVars_, bcTypes_); - dest -= localResidual().residual(); - destStorage -= localResidual().storageTerm()[scvIdx]; + partialDeriv -= localResidual().residual(); + storageDeriv -= localResidual().storageTerm()[scvIdx]; } else { // we are using forward differences, i.e. we don't need to // calculate f(x - \epsilon) and we can recycle the // (already calculated) residual f(x) - dest -= residual_; - destStorage -= storageTerm_[scvIdx]; + partialDeriv -= residual_; + storageDeriv -= storageTerm_[scvIdx]; } // divide difference in residuals by the magnitude of the // deflections between the two function evaluation - dest /= delta; - destStorage /= delta; + partialDeriv /= delta; + storageDeriv /= delta; // restore the original state of the element's volume variables curVolVars_[scvIdx] = origVolVars; #if HAVE_VALGRIND - for (unsigned i = 0; i < dest.size(); ++i) - Valgrind::CheckDefined(dest[i]); + for (unsigned i = 0; i < partialDeriv.size(); ++i) + Valgrind::CheckDefined(partialDeriv[i]); #endif } @@ -486,7 +496,7 @@ protected: */ void updateLocalJacobian_(int scvIdx, int pvIdx, - const ElementSolutionVector &deriv, + const ElementSolutionVector &partialDeriv, const PrimaryVariables &storageDeriv) { // store the derivative of the storage term @@ -496,7 +506,7 @@ protected: for (int i = 0; i < fvElemGeom_.numVertices; i++) { - if (jacAsm_().vertexColor(elem_(), i) == Green) { + if (jacAsm_().vertexColor(element_(), i) == Green) { // Green vertices are not to be changed! continue; } @@ -506,7 +516,7 @@ protected: // the residual of equation 'eqIdx' at vertex 'i' // depending on the primary variable 'pvIdx' at vertex // 'scvIdx'. - this->A_[i][scvIdx][eqIdx][pvIdx] = deriv[i][eqIdx]; + this->A_[i][scvIdx][eqIdx][pvIdx] = partialDeriv[i][eqIdx]; Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]); } } diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh index 5e296ab3a31bac53f11f24377c313fca85442912..d4b8531617899a6a0b0c355178420f3959f71771 100644 --- a/dumux/boxmodels/common/boxlocalresidual.hh +++ b/dumux/boxmodels/common/boxlocalresidual.hh @@ -92,11 +92,11 @@ public: * This assumes that all objects of the simulation have been fully * allocated but not necessarily initialized completely. * - * \param prob The representation of the physical problem to be + * \param problem The representation of the physical problem to be * solved. */ - void init(Problem &prob) - { problemPtr_ = &prob; } + void init(Problem &problem) + { problemPtr_ = &problem; } /*! * \brief Compute the local residual, i.e. the deviation of the @@ -107,10 +107,10 @@ public: */ void eval(const Element &element) { - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeometry; - fvElemGeom.update(gridView_(), element); - fvElemGeomPtr_ = &fvElemGeom; + fvGeometry.update(gridView_(), element); + fvElemGeomPtr_ = &fvGeometry; ElementVolumeVariables volVarsPrev, volVarsCur; // update the hints @@ -151,9 +151,9 @@ public: { elemPtr_ = &element; - FVElementGeometry fvElemGeom; - fvElemGeom.update(gridView_(), element); - fvElemGeomPtr_ = &fvElemGeom; + FVElementGeometry fvGeometry; + fvGeometry.update(gridView_(), element); + fvElemGeomPtr_ = &fvGeometry; ElementBoundaryTypes bcTypes; bcTypes.update(problem_(), element, fvGeometry_()); @@ -187,9 +187,9 @@ public: { elemPtr_ = &element; - FVElementGeometry fvElemGeom; - fvElemGeom.update(gridView_(), element); - fvElemGeomPtr_ = &fvElemGeom; + FVElementGeometry fvGeometry; + fvGeometry.update(gridView_(), element); + fvElemGeomPtr_ = &fvGeometry; ElementBoundaryTypes bcTypes; bcTypes.update(problem_(), element, fvGeometry_()); @@ -209,7 +209,7 @@ public: * * \param element The DUNE Codim<0> entity for which the residual * ought to be calculated - * \param fvElemGeom The finite-volume geometry of the element + * \param fvGeometry The finite-volume geometry of the element * \param prevVolVars The volume averaged variables for all * sub-control volumes of the element at the previous * time level @@ -220,7 +220,7 @@ public: * vertices of the element */ void eval(const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, const ElementVolumeVariables &prevVolVars, const ElementVolumeVariables &curVolVars, const ElementBoundaryTypes &bcTypes) @@ -229,14 +229,14 @@ public: Valgrind::CheckDefined(curVolVars); #if !defined NDEBUG && HAVE_VALGRIND - for (int i=0; i < fvElemGeom.numVertices; i++) { + for (int i=0; i < fvGeometry.numVertices; i++) { prevVolVars[i].checkDefined(); curVolVars[i].checkDefined(); } #endif // HAVE_VALGRIND elemPtr_ = &element; - fvElemGeomPtr_ = &fvElemGeom; + fvElemGeomPtr_ = &fvGeometry; bcTypesPtr_ = &bcTypes; prevVolVarsPtr_ = &prevVolVars; curVolVarsPtr_ = &curVolVars; @@ -353,16 +353,16 @@ protected: */ void evalDirichlet_() { - PrimaryVariables tmp(0); - for (int i = 0; i < fvGeometry_().numVertices; ++i) { - const BoundaryTypes &bcTypes = bcTypes_(i); + PrimaryVariables dirichletValues(0); + for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) { + const BoundaryTypes &bcTypes = bcTypes_(scvIdx); if (! bcTypes.hasDirichlet()) continue; // ask the problem for the dirichlet values - const VertexPointer vPtr = element_().template subEntity<dim>(i); - Valgrind::SetUndefined(tmp); - asImp_().problem_().dirichlet(tmp, *vPtr); + const VertexPointer vPtr = element_().template subEntity<dim>(scvIdx); + Valgrind::SetUndefined(dirichletValues); + asImp_().problem_().dirichlet(dirichletValues, *vPtr); // set the dirichlet conditions for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { @@ -370,12 +370,12 @@ protected: continue; int pvIdx = bcTypes.eqToDirichletIndex(eqIdx); assert(0 <= pvIdx && pvIdx < numEq); - Valgrind::CheckDefined(tmp[pvIdx]); + Valgrind::CheckDefined(dirichletValues[pvIdx]); - residual_[i][eqIdx] = - curPrimaryVar_(i, pvIdx) - tmp[pvIdx]; + residual_[scvIdx][eqIdx] = + curPrimaryVar_(scvIdx, pvIdx) - dirichletValues[pvIdx]; - storageTerm_[i][eqIdx] = 0.0; + storageTerm_[scvIdx][eqIdx] = 0.0; }; }; } @@ -405,7 +405,7 @@ protected: faceVertIdx < numFaceVerts; ++faceVertIdx) { - int elemVertIdx = refElem.subEntity(faceIdx, + int scvIdx = refElem.subEntity(faceIdx, 1, faceVertIdx, dim); @@ -416,13 +416,13 @@ protected: // add the residual of all vertices of the boundary // segment asImp_().evalNeumannSegment_(isIt, - elemVertIdx, + scvIdx, boundaryFaceIdx); // evaluate the outflow conditions at the boundary face // ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model // available and not thoroughly tested. asImp_().evalOutflowSegment(isIt, - elemVertIdx, + scvIdx, boundaryFaceIdx); } } @@ -437,29 +437,29 @@ protected: int boundaryFaceIdx) { // temporary vector to store the neumann boundary fluxes - PrimaryVariables values(0.0); + PrimaryVariables neumannFlux(0.0); const BoundaryTypes &bcTypes = bcTypes_(scvIdx); // deal with neumann boundaries if (bcTypes.hasNeumann()) { - Valgrind::SetUndefined(values); - problem_().boxSDNeumann(values, + Valgrind::SetUndefined(neumannFlux); + problem_().boxSDNeumann(neumannFlux, element_(), fvGeometry_(), *isIt, scvIdx, boundaryFaceIdx, curVolVars_()); - values *= + neumannFlux *= fvGeometry_().boundaryFace[boundaryFaceIdx].area * curVolVars_(scvIdx).extrusionFactor(); - Valgrind::CheckDefined(values); + Valgrind::CheckDefined(neumannFlux); // set the neumann conditions for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { if (!bcTypes.isNeumann(eqIdx)) continue; - residual_[scvIdx][eqIdx] += values[eqIdx]; + residual_[scvIdx][eqIdx] += neumannFlux[eqIdx]; } } } @@ -472,15 +472,15 @@ protected: { // calculate the mass flux over the faces and subtract // it from the local rates - for (int k = 0; k < fvGeometry_().numEdges; k++) + for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++) { - int i = fvGeometry_().subContVolFace[k].i; - int j = fvGeometry_().subContVolFace[k].j; + int i = fvGeometry_().subContVolFace[scvfIdx].i; + int j = fvGeometry_().subContVolFace[scvfIdx].j; PrimaryVariables flux; Valgrind::SetUndefined(flux); - asImp_().computeFlux(flux, k); + asImp_().computeFlux(flux, scvfIdx); Valgrind::CheckDefined(flux); Scalar extrusionFactor = @@ -519,13 +519,13 @@ protected: // calculate the amount of conservation each quantity inside // all sub control volumes - for (int i=0; i < fvGeometry_().numVertices; i++) { - Valgrind::SetUndefined(storageTerm_[i]); - asImp_().computeStorage(storageTerm_[i], i, /*isOldSol=*/false); - storageTerm_[i] *= - fvGeometry_().subContVol[i].volume - * curVolVars_(i).extrusionFactor(); - Valgrind::CheckDefined(storageTerm_[i]); + for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) { + Valgrind::SetUndefined(storageTerm_[scvIdx]); + asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, /*isOldSol=*/false); + storageTerm_[scvIdx] *= + fvGeometry_().subContVol[scvIdx].volume + * curVolVars_(scvIdx).extrusionFactor(); + Valgrind::CheckDefined(storageTerm_[scvIdx]); } } @@ -537,12 +537,12 @@ protected: void evalVolumeTerms_() { // evaluate the volume terms (storage + source terms) - for (int i=0; i < fvGeometry_().numVertices; i++) + for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) { Scalar extrusionFactor = - curVolVars_(i).extrusionFactor(); + curVolVars_(scvIdx).extrusionFactor(); - PrimaryVariables tmp(0.); + PrimaryVariables values(0.0); // mass balance within the element. this is the // \f$\frac{m}{\partial t}\f$ term if using implicit @@ -550,30 +550,30 @@ protected: // // TODO (?): we might need a more explicit way for // doing the time discretization... - Valgrind::SetUndefined(storageTerm_[i]); - Valgrind::SetUndefined(tmp); - asImp_().computeStorage(storageTerm_[i], i, false); - asImp_().computeStorage(tmp, i, true); - Valgrind::CheckDefined(storageTerm_[i]); - Valgrind::CheckDefined(tmp); - - storageTerm_[i] -= tmp; - storageTerm_[i] *= - fvGeometry_().subContVol[i].volume + Valgrind::SetUndefined(storageTerm_[scvIdx]); + Valgrind::SetUndefined(values); + asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, false); + asImp_().computeStorage(values, scvIdx, true); + Valgrind::CheckDefined(storageTerm_[scvIdx]); + Valgrind::CheckDefined(values); + + storageTerm_[scvIdx] -= values; + storageTerm_[scvIdx] *= + fvGeometry_().subContVol[scvIdx].volume / problem_().timeManager().timeStepSize() * extrusionFactor; - residual_[i] += storageTerm_[i]; + residual_[scvIdx] += storageTerm_[scvIdx]; // subtract the source term from the local rate - Valgrind::SetUndefined(tmp); - asImp_().computeSource(tmp, i); - Valgrind::CheckDefined(tmp); - tmp *= fvGeometry_().subContVol[i].volume * extrusionFactor; - residual_[i] -= tmp; + Valgrind::SetUndefined(values); + asImp_().computeSource(values, scvIdx); + Valgrind::CheckDefined(values); + values *= fvGeometry_().subContVol[scvIdx].volume * extrusionFactor; + residual_[scvIdx] -= values; // make sure that only defined quantities were used // to calculate the residual. - Valgrind::CheckDefined(residual_[i]); + Valgrind::CheckDefined(residual_[scvIdx]); } } diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh index dd0e33b3f2f4779e349893d914ead592df148f71..d7c0c6c16fac725a0b62a3696ef20affb5a125fb 100644 --- a/dumux/boxmodels/common/boxmodel.hh +++ b/dumux/boxmodels/common/boxmodel.hh @@ -103,12 +103,12 @@ public: /*! * \brief Apply the initial conditions to the model. * - * \param prob The object representing the problem which needs to + * \param problem The object representing the problem which needs to * be simulated. */ - void init(Problem &prob) + void init(Problem &problem) { - problemPtr_ = &prob; + problemPtr_ = &problem; updateBoundaryIndices_(); @@ -139,18 +139,18 @@ public: uPrev_ = uCur_; } - void setHints(const Element &elem, + void setHints(const Element &element, ElementVolumeVariables &prevVolVars, ElementVolumeVariables &curVolVars) const { if (!enableHints_) return; - int n = elem.template count<dim>(); + int n = element.template count<dim>(); prevVolVars.resize(n); curVolVars.resize(n); for (int i = 0; i < n; ++i) { - int globalIdx = vertexMapper().map(elem, i, dim); + int globalIdx = vertexMapper().map(element, i, dim); if (!hintsUsable_[globalIdx]) { curVolVars[i].setHint(NULL); @@ -163,16 +163,16 @@ public: } }; - void setHints(const Element &elem, + void setHints(const Element &element, ElementVolumeVariables &curVolVars) const { if (!enableHints_) return; - int n = elem.template count<dim>(); + int n = element.template count<dim>(); curVolVars.resize(n); for (int i = 0; i < n; ++i) { - int globalIdx = vertexMapper().map(elem, i, dim); + int globalIdx = vertexMapper().map(element, i, dim); if (!hintsUsable_[globalIdx]) curVolVars[i].setHint(NULL); @@ -189,17 +189,17 @@ public: prevHints_ = curHints_; }; - void updateCurHints(const Element &elem, - const ElementVolumeVariables &ev) const + void updateCurHints(const Element &element, + const ElementVolumeVariables &elemVolVars) const { if (!enableHints_) return; - for (int i = 0; i < ev.size(); ++i) { - int globalIdx = vertexMapper().map(elem, i, dim); - curHints_[globalIdx] = ev[i]; + for (int i = 0; i < elemVolVars.size(); ++i) { + int globalIdx = vertexMapper().map(element, i, dim); + curHints_[globalIdx] = elemVolVars[i]; if (!hintsUsable_[globalIdx]) - prevHints_[globalIdx] = ev[i]; + prevHints_[globalIdx] = elemVolVars[i]; hintsUsable_[globalIdx] = true; } }; @@ -209,15 +209,15 @@ public: * \brief Compute the global residual for an arbitrary solution * vector. * - * \param dest Stores the result + * \param residual Stores the result * \param u The solution for which the residual ought to be calculated */ - Scalar globalResidual(SolutionVector &dest, + Scalar globalResidual(SolutionVector &residual, const SolutionVector &u) { SolutionVector tmp(curSol()); curSol() = u; - Scalar res = globalResidual(dest); + Scalar res = globalResidual(residual); curSol() = tmp; return res; } @@ -226,11 +226,11 @@ public: * \brief Compute the global residual for the current solution * vector. * - * \param dest Stores the result + * \param residual Stores the result */ - Scalar globalResidual(SolutionVector &dest) + Scalar globalResidual(SolutionVector &residual) { - dest = 0; + residual = 0; ElementIterator elemIt = gridView_().template begin<0>(); const ElementIterator elemEndIt = gridView_().template end<0>(); @@ -239,19 +239,19 @@ public: for (int i = 0; i < elemIt->template count<dim>(); ++i) { int globalI = vertexMapper().map(*elemIt, i, dim); - dest[globalI] += localResidual().residual(i); + residual[globalI] += localResidual().residual(i); } }; // calculate the square norm of the residual - Scalar result2 = dest.two_norm2(); + Scalar result2 = residual.two_norm2(); if (gridView_().comm().size() > 1) result2 = gridView_().comm().sum(result2); // add up the residuals on the process borders if (gridView_().comm().size() > 1) { VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper> - sumHandle(dest, vertexMapper()); + sumHandle(residual, vertexMapper()); gridView_().communicate(sumHandle, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication); @@ -264,11 +264,11 @@ public: * \brief Compute the integral over the domain of the storage * terms of all conservation quantities. * - * \param dest Stores the result + * \param storage Stores the result */ - void globalStorage(PrimaryVariables &dest) + void globalStorage(PrimaryVariables &storage) { - dest = 0; + storage = 0; ElementIterator elemIt = gridView_().template begin<0>(); const ElementIterator elemEndIt = gridView_().template end<0>(); @@ -276,11 +276,11 @@ public: localResidual().evalStorage(*elemIt); for (int i = 0; i < elemIt->template count<dim>(); ++i) - dest += localResidual().storageTerm()[i]; + storage += localResidual().storageTerm()[i]; }; if (gridView_().comm().size() > 1) - dest = gridView_().comm().sum(dest); + storage = gridView_().comm().sum(storage); } /*! @@ -513,13 +513,13 @@ public: * * \param outstream The stream into which the vertex data should * be serialized to - * \param vert The DUNE Codim<dim> entity which's data should be + * \param vertex The DUNE Codim<dim> entity which's data should be * serialized */ void serializeEntity(std::ostream &outstream, - const Vertex &vert) + const Vertex &vertex) { - int vertIdx = dofMapper().map(vert); + int vertIdx = dofMapper().map(vertex); // write phase state if (!outstream.good()) { @@ -539,13 +539,13 @@ public: * * \param instream The stream from which the vertex data should * be deserialized from - * \param vert The DUNE Codim<dim> entity which's data should be + * \param vertex The DUNE Codim<dim> entity which's data should be * deserialized */ void deserializeEntity(std::istream &instream, - const Vertex &vert) + const Vertex &vertex) { - int vertIdx = dofMapper().map(vert); + int vertIdx = dofMapper().map(vertex); for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { if (!instream.good()) DUNE_THROW(Dune::IOError, @@ -621,21 +621,20 @@ public: { typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField; - SolutionVector globalResid(u); - asImp_().globalResidual(globalResid, u); + SolutionVector residual(u); + asImp_().globalResidual(residual, u); // create the required scalar fields unsigned numVertices = this->gridView_().size(dim); - //unsigned numElements = this->gridView_().size(0); // global defect of the two auxiliary equations ScalarField* def[numEq]; ScalarField* delta[numEq]; ScalarField* x[numEq]; - for (int i = 0; i < numEq; ++i) { - x[i] = writer.allocateManagedBuffer(numVertices); - delta[i] = writer.allocateManagedBuffer(numVertices); - def[i] = writer.allocateManagedBuffer(numVertices); + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + x[eqIdx] = writer.allocateManagedBuffer(numVertices); + delta[eqIdx] = writer.allocateManagedBuffer(numVertices); + def[eqIdx] = writer.allocateManagedBuffer(numVertices); } VertexIterator vIt = this->gridView_().template begin<dim>(); @@ -643,21 +642,21 @@ public: for (; vIt != vEndIt; ++ vIt) { int globalIdx = vertexMapper().map(*vIt); - for (int i = 0; i < numEq; ++i) { - (*x[i])[globalIdx] = u[globalIdx][i]; - (*delta[i])[globalIdx] = - deltaU[globalIdx][i]; - (*def[i])[globalIdx] = globalResid[globalIdx][i]; + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + (*x[eqIdx])[globalIdx] = u[globalIdx][eqIdx]; + (*delta[eqIdx])[globalIdx] = - deltaU[globalIdx][eqIdx]; + (*def[eqIdx])[globalIdx] = residual[globalIdx][eqIdx]; } } - for (int i = 0; i < numEq; ++i) { + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { std::ostringstream oss; - oss.str(""); oss << "x_" << i; - writer.attachVertexData(*x[i], oss.str()); - oss.str(""); oss << "delta_" << i; - writer.attachVertexData(*delta[i], oss.str()); - oss.str(""); oss << "defect_" << i; - writer.attachVertexData(*def[i], oss.str()); + oss.str(""); oss << "x_" << eqIdx; + writer.attachVertexData(*x[eqIdx], oss.str()); + oss.str(""); oss << "delta_" << eqIdx; + writer.attachVertexData(*delta[eqIdx], oss.str()); + oss.str(""); oss << "defect_" << eqIdx; + writer.attachVertexData(*def[eqIdx], oss.str()); } asImp_().addOutputVtkFields(u, writer); @@ -686,8 +685,8 @@ public: // global defect of the two auxiliary equations ScalarField* x[numEq]; - for (int i = 0; i < numEq; ++i) { - x[i] = writer.allocateManagedBuffer(numVertices); + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + x[eqIdx] = writer.allocateManagedBuffer(numVertices); } VertexIterator vIt = this->gridView_().template begin<dim>(); @@ -695,15 +694,15 @@ public: for (; vIt != vEndIt; ++ vIt) { int globalIdx = vertexMapper().map(*vIt); - for (int i = 0; i < numEq; ++i) { - (*x[i])[globalIdx] = sol[globalIdx][i]; + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { + (*x[eqIdx])[globalIdx] = sol[globalIdx][eqIdx]; } } - for (int i = 0; i < numEq; ++i) { + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) { std::ostringstream oss; - oss << "primaryVar_" << i; - writer.attachVertexData(*x[i], oss.str()); + oss << "primaryVar_" << eqIdx; + writer.attachVertexData(*x[eqIdx], oss.str()); } } @@ -727,12 +726,12 @@ public: * \brief Returns true if a vertex is located on the grid's * boundary. * - * \param elem A DUNE Codim<0> entity which contains the control + * \param element A DUNE Codim<0> entity which contains the control * volume's associated vertex. - * \param vIdx The local vertex index inside elem + * \param vIdx The local vertex index inside element */ - bool onBoundary(const Element &elem, int vIdx) const - { return onBoundary(vertexMapper().map(elem, vIdx, dim)); } + bool onBoundary(const Element &element, int vIdx) const + { return onBoundary(vertexMapper().map(element, vIdx, dim)); } /*! * \brief Fill the fluid state according to the primary variables. @@ -792,24 +791,24 @@ protected: uCur_ = Scalar(0.0); boxVolume_ = Scalar(0.0); - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeometry; // iterate through leaf grid and evaluate initial // condition at the center of each sub control volume // // TODO: the initial condition needs to be unique for // each vertex. we should think about the API... - ElementIterator it = gridView_().template begin<0>(); - const ElementIterator &eendit = gridView_().template end<0>(); - for (; it != eendit; ++it) { + ElementIterator eIt = gridView_().template begin<0>(); + const ElementIterator &eEndIt = gridView_().template end<0>(); + for (; eIt != eEndIt; ++eIt) { // deal with the current element - fvElemGeom.update(gridView_(), *it); + fvGeometry.update(gridView_(), *eIt); // loop over all element vertices, i.e. sub control volumes - for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; scvIdx++) + for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; scvIdx++) { // map the local vertex index to the global one - int globalIdx = vertexMapper().map(*it, + int globalIdx = vertexMapper().map(*eIt, scvIdx, dim); @@ -818,8 +817,8 @@ protected: PrimaryVariables initVal; Valgrind::SetUndefined(initVal); problem_().initial(initVal, - *it, - fvElemGeom, + *eIt, + fvGeometry, scvIdx); Valgrind::CheckDefined(initVal); @@ -827,8 +826,8 @@ protected: // volumes. If the initial values disagree for // different sub control volumes, the initial value // will be the arithmetic mean. - initVal *= fvElemGeom.subContVol[scvIdx].volume; - boxVolume_[globalIdx] += fvElemGeom.subContVol[scvIdx].volume; + initVal *= fvGeometry.subContVol[scvIdx].volume; + boxVolume_[globalIdx] += fvGeometry.subContVol[scvIdx].volume; uCur_[globalIdx] += initVal; Valgrind::CheckDefined(uCur_[globalIdx]); } diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh index a9665948b77f9d06baacea271d6ec2b0e9bb2685..cab2c1873bba554ea3af79dd5ccdc9520b35a67c 100644 --- a/dumux/boxmodels/common/boxproblem.hh +++ b/dumux/boxmodels/common/boxproblem.hh @@ -222,7 +222,7 @@ public: * * \param values The neumann values for the conservation equations [kg / (m^2 *s )] * \param element The finite element - * \param fvElemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param is The intersection between element and boundary * \param scvIdx The local vertex index * \param boundaryFaceIdx The index of the boundary face @@ -233,7 +233,7 @@ public: */ void boxSDNeumann(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, const Intersection &is, int scvIdx, int boundaryFaceIdx, @@ -242,7 +242,7 @@ public: // forward it to the interface without the volume variables asImp_().neumann(values, element, - fvElemGeom, + fvGeometry, is, scvIdx, boundaryFaceIdx); @@ -254,7 +254,7 @@ public: * * \param values The neumann values for the conservation equations [kg / (m^2 *s )] * \param element The finite element - * \param fvElemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param is The intersection between element and boundary * \param scvIdx The local vertex index * \param boundaryFaceIdx The index of the boundary face @@ -264,13 +264,13 @@ public: */ void neumann(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, const Intersection &is, int scvIdx, int boundaryFaceIdx) const { // forward it to the interface with only the global position - asImp_().neumannAtPos(values, fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal); + asImp_().neumannAtPos(values, fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal); } /*! @@ -304,7 +304,7 @@ public: * * \param values The source and sink values for the conservation equations * \param element The finite element - * \param fvElemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param scvIdx The local vertex index * \param elemVolVars All volume variables for the element * @@ -314,12 +314,12 @@ public: */ void boxSDSource(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, int scvIdx, const ElementVolumeVariables &elemVolVars) const { // forward to solution independent, box specific interface - asImp_().source(values, element, fvElemGeom, scvIdx); + asImp_().source(values, element, fvGeometry, scvIdx); } /*! @@ -328,7 +328,7 @@ public: * * \param values The source and sink values for the conservation equations * \param element The finite element - * \param fvElemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param scvIdx The local vertex index * * For this method, the \a values parameter stores the rate mass @@ -337,11 +337,11 @@ public: */ void source(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, int scvIdx) const { // forward to generic interface - asImp_().sourceAtPos(values, fvElemGeom.subContVol[scvIdx].global); + asImp_().sourceAtPos(values, fvGeometry.subContVol[scvIdx].global); } /*! @@ -370,7 +370,7 @@ public: * * \param values The initial values for the primary variables * \param element The finite element - * \param fvElemGeom The finite-volume geometry in the box scheme + * \param fvGeometry The finite-volume geometry in the box scheme * \param scvIdx The local vertex index * * For this method, the \a values parameter stores primary @@ -378,11 +378,11 @@ public: */ void initial(PrimaryVariables &values, const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, int scvIdx) const { // forward to generic interface - asImp_().initialAtPos(values, fvElemGeom.subContVol[scvIdx].global); + asImp_().initialAtPos(values, fvGeometry.subContVol[scvIdx].global); } /*! @@ -415,11 +415,11 @@ public: * are assumed to extend 1 m to the back. */ Scalar boxExtrusionFactor(const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, int scvIdx) const { // forward to generic interface - return asImp_().extrusionFactorAtPos(fvElemGeom.subContVol[scvIdx].global); + return asImp_().extrusionFactorAtPos(fvGeometry.subContVol[scvIdx].global); } /*! diff --git a/dumux/boxmodels/common/boxvolumevariables.hh b/dumux/boxmodels/common/boxvolumevariables.hh index da479af3f272dc9737716ce33a1ed26692a25b5b..f5e0d60ff3f6f09e675a9eb0f062e9b97dfb1050 100644 --- a/dumux/boxmodels/common/boxvolumevariables.hh +++ b/dumux/boxmodels/common/boxvolumevariables.hh @@ -112,7 +112,7 @@ public: * \param problem The object specifying the problem which ought to * be simulated * \param element An element which contains part of the control volume - * \param elemGeom The finite volume geometry for the element + * \param fvGeometry The finite volume geometry for the element * \param scvIdx Local index of the sub control volume which is inside the element * \param isOldSol Specifies whether this is the previous solution or the current one * @@ -124,13 +124,13 @@ public: void update(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, + const FVElementGeometry &fvGeometry, int scvIdx, bool isOldSol) { Valgrind::CheckDefined(priVars); primaryVars_ = priVars; - extrusionFactor_ = problem.boxExtrusionFactor(element, elemGeom, scvIdx); + extrusionFactor_ = problem.boxExtrusionFactor(element, fvGeometry, scvIdx); } /*! diff --git a/dumux/boxmodels/common/porousmediaboxproblem.hh b/dumux/boxmodels/common/porousmediaboxproblem.hh index 2367e47e3af8d4308c7afbfa7aa3118a76a3df8c..e1390a595271b78740c164125e8535461102e89c 100644 --- a/dumux/boxmodels/common/porousmediaboxproblem.hh +++ b/dumux/boxmodels/common/porousmediaboxproblem.hh @@ -63,7 +63,7 @@ class PorousMediaBoxProblem : public BoxProblem<TypeTag> typedef typename GridView::ctype CoordScalar; typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef Dune::FieldVector<Scalar, dim> Vector; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: /*! @@ -105,13 +105,13 @@ public: * * \param element The DUNE Codim<0> enitiy which intersects with * the finite volume. - * \param fvGeom The finite volume geometry of the element. + * \param fvGeometry The finite volume geometry of the element. * \param scvIdx The local index of the sub control volume inside the element */ Scalar boxTemperature(const Element &element, - const FVElementGeometry fvGeom, + const FVElementGeometry fvGeometry, int scvIdx) const - { return asImp_().temperatureAtPos(fvGeom.subContVol[scvIdx].global); } + { return asImp_().temperatureAtPos(fvGeometry.subContVol[scvIdx].global); } /*! * \brief Returns the temperature \f$\mathrm{[K]}\f$ at a given global position. @@ -140,10 +140,10 @@ public: * This is the box discretization specific interface. By default * it just calls gravityAtPos(). */ - const Vector &boxGravity(const Element &element, - const FVElementGeometry &fvGeom, + const DimVector &boxGravity(const Element &element, + const FVElementGeometry &fvGeometry, int scvIdx) const - { return asImp_().gravityAtPos(fvGeom.subContVol[scvIdx].global); } + { return asImp_().gravityAtPos(fvGeometry.subContVol[scvIdx].global); } /*! * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. @@ -151,7 +151,7 @@ public: * This is discretization independent interface. By default it * just calls gravity(). */ - const Vector &gravityAtPos(const GlobalPosition &pos) const + const DimVector &gravityAtPos(const GlobalPosition &pos) const { return asImp_().gravity(); } /*! @@ -163,7 +163,7 @@ public: * property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds, * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. */ - const Vector &gravity() const + const DimVector &gravity() const { return gravity_; } /*! @@ -196,7 +196,7 @@ protected: const Implementation &asImp_() const { return *static_cast<const Implementation *>(this); } - Vector gravity_; + DimVector gravity_; // fluids and material properties SpatialParams* spatialParams_;