diff --git a/configure.ac b/configure.ac
index 1445dfef0bf5f03200ff8cd0640861b8a956083f..85f211277d43fdc2e10f7e3b7474385c8a5a3019 100644
--- a/configure.ac
+++ b/configure.ac
@@ -82,8 +82,9 @@ AC_CONFIG_FILES([dumux.pc
dumux/material/constraintsolvers/Makefile
dumux/material/eos/Makefile
dumux/multidomain/Makefile
- dumux/multidomain/common/Makefile
dumux/multidomain/2cstokes2p2c/Makefile
+ dumux/multidomain/common/Makefile
+ dumux/multidomain/couplinglocalresiduals/Makefile
dumux/nonlinear/Makefile
dumux/parallel/Makefile
m4/Makefile
diff --git a/dumux/multidomain/Makefile.am b/dumux/multidomain/Makefile.am
index 4b6c342f0c47e220af49cb3ddf73c66a44eda26c..2f533181d02c88fc432bed6c9fc4c2f402de058b 100644
--- a/dumux/multidomain/Makefile.am
+++ b/dumux/multidomain/Makefile.am
@@ -1,4 +1,6 @@
-SUBDIRS = common \
+SUBDIRS = \
+ common \
+ couplinglocalresiduals \
2cstokes2p2c
multidomaindir = $(includedir)/dumux/multidomain
diff --git a/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7c10fcb12c2cc078a9961992225072800b5b5e12
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
@@ -0,0 +1,343 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Supplements the TwoPTwoCLocalResidual by the required functions for a coupled application.
+ *
+ */
+#ifndef DUMUX_2P2C_COUPLING_LOCAL_RESIDUAL_HH
+#define DUMUX_2P2C_COUPLING_LOCAL_RESIDUAL_HH
+
+#include <dumux/implicit/2p2c/2p2clocalresidual.hh>
+
+namespace Dumux
+{
+
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Supplements the TwoPTwoCModel by the required functions for a coupled application.
+ */
+template<class TypeTag>
+class TwoPTwoCCouplingLocalResidual : public TwoPTwoCLocalResidual<TypeTag>
+{
+ typedef TwoPTwoCCouplingLocalResidual<TypeTag> ThisType;
+ typedef TwoPTwoCLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ enum { dim = GridView::dimension };
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum {
+ pressureIdx = Indices::pressureIdx
+ };
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+ enum {
+ massBalanceIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
+ contiWEqIdx = Indices::contiWEqIdx
+ };
+ enum {
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx
+ };
+ enum { phaseIdx = nPhaseIdx }; // index of the phase for the phase flux calculation
+ enum { compIdx = wCompIdx}; // index of the component for the phase flux calculation
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar,1> > ElementFluxVector;
+
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ TwoPTwoCCouplingLocalResidual()
+ { };
+
+ void evalBoundary_()
+ {
+ ParentType::evalBoundary_();
+
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+ const ReferenceElement &refElement = ReferenceElements::general(this->element_().geometry().type());
+
+ // evaluate Dirichlet-like coupling conditions
+ for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
+ {
+ // evaluate boundary conditions for the intersections of
+ // the current element
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
+ for (; isIt != endIt; ++isIt)
+ {
+ // handle only intersections on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ // assemble the boundary for all vertices of the current face
+ const int faceIdx = isIt->indexInInside();
+ const int numFaceVertices = refElement.size(faceIdx, 1, dim);
+
+ // loop over the single vertices on the current face
+ for (int faceVertIdx = 0; faceVertIdx < numFaceVertices; ++faceVertIdx)
+ {
+ const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ const int elemVertIdx = refElement.subEntity(faceIdx, 1, faceVertIdx, dim);
+ // only evaluate, if we consider the same face vertex as in the outer
+ // loop over the element vertices
+ if (elemVertIdx != idx)
+ continue;
+
+ //for the corner points, the boundary flux across the vertical non-coupling boundary faces
+ //has to be calculated to fulfill the mass balance
+ //convert suddomain intersection into multidomain intersection and check whether it is an outer boundary
+ if(!GridView::Grid::multiDomainIntersection(*isIt).neighbor()
+ && this->boundaryHasMortarCoupling_(this->bcTypes_(idx)))//this->boundaryHasNeumann_(this->bcTypes_(idx)))
+ {
+ const DimVector& globalPos = this->fvGeometry_().subContVol[idx].global;
+ //problem specific function, in problem orientation of interface is known
+ if(this->problem_().isInterfaceCornerPoint(globalPos))
+ {
+ //check whether the second boundary node on the vertical edge is Neumann-Null
+ const int numVertices = refElement.size(dim);
+ bool evalBoundaryFlux = false;
+ for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
+ {
+ for(int i= 0; i < numVertices; i++)
+ {
+ //if vertex is on boundary and not the coupling vertex: check whether an outflow condition is set
+ if(this->model_().onBoundary(this->element_(), i) && i!=elemVertIdx)
+ if (!this->bcTypes_(i).isOutflow(equationIdx))
+ evalBoundaryFlux = true;
+ }
+
+ PrimaryVariables values(0.0);
+ //calculate the actual boundary fluxes and add to residual
+ this->asImp_()->computeFlux(values, boundaryFaceIdx, true /*on boundary*/);
+ if(evalBoundaryFlux)
+ this->residual_[idx][equationIdx] += values[equationIdx];
+ }
+
+ }
+ }
+
+ if (boundaryHasCoupling_(this->bcTypes_(idx)))
+ evalCouplingVertex_(idx);
+ }
+ }
+ }
+ }
+
+
+ Scalar evalPhaseStorage(const int scvIdx) const
+ {
+ Scalar phaseStorage = computePhaseStorage(scvIdx, false);
+ Scalar oldPhaseStorage = computePhaseStorage(scvIdx, true);;
+ Valgrind::CheckDefined(phaseStorage);
+ Valgrind::CheckDefined(oldPhaseStorage);
+
+ phaseStorage -= oldPhaseStorage;
+ phaseStorage *= this->fvGeometry_().subContVol[scvIdx].volume
+ / this->problem_().timeManager().timeStepSize()
+ * this->curVolVars_(scvIdx).extrusionFactor();
+ Valgrind::CheckDefined(phaseStorage);
+
+ return phaseStorage;
+ }
+
+ /*!
+ * compute storage term of all components within all phases
+ */
+ Scalar computePhaseStorage(const int scvIdx, bool usePrevSol) const
+ {
+ const ElementVolumeVariables &elemVolVars =
+ usePrevSol ? this->prevVolVars_() : this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+ return volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFraction(phaseIdx, compIdx)
+ * volVars.porosity();
+ }
+
+ /*!
+ * \brief Compute the fluxes within the different fluid phases. This is
+ * merely for the computation of flux output.
+ */
+ void evalPhaseFluxes()
+ {
+ elementFluxes_.resize(this->fvGeometry_().numScv);
+ elementFluxes_ = 0.;
+
+ Scalar flux(0.);
+
+ // calculate the mass flux over the faces and subtract
+ // it from the local rates
+ for (int faceIdx = 0; faceIdx < this->fvGeometry_().numScvf; faceIdx++)
+ {
+ FluxVariables vars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ faceIdx,
+ this->curVolVars_());
+
+ int i = this->fvGeometry_().subContVolFace[faceIdx].i;
+ int j = this->fvGeometry_().subContVolFace[faceIdx].j;
+
+ const Scalar extrusionFactor =
+ (this->curVolVars_(i).extrusionFactor()
+ + this->curVolVars_(j).extrusionFactor())
+ / 2;
+ flux = computeAdvectivePhaseFlux(vars) +
+ computeDiffusivePhaseFlux(vars);
+ flux *= extrusionFactor;
+
+ elementFluxes_[i] += flux;
+ elementFluxes_[j] -= flux;
+ }
+ }
+
+ /*!
+ * \brief Compute the advective fluxes within the different phases.
+ */
+ Scalar computeAdvectivePhaseFlux(const FluxVariables &fluxVars) const
+ {
+ Scalar advectivePhaseFlux = 0.;
+ const Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up =
+ this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn =
+ this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
+ if (massUpwindWeight > 0.0)
+ // upstream vertex
+ advectivePhaseFlux +=
+ fluxVars.volumeFlux(phaseIdx)
+ * massUpwindWeight
+ * up.density(phaseIdx)
+ * up.fluidState().massFraction(phaseIdx, compIdx);
+ if (massUpwindWeight < 1.0)
+ // downstream vertex
+ advectivePhaseFlux +=
+ fluxVars.volumeFlux(phaseIdx)
+ * (1 - massUpwindWeight)
+ * dn.density(phaseIdx)
+ * dn.fluidState().massFraction(phaseIdx, compIdx);
+
+ return advectivePhaseFlux;
+ }
+
+ /*!
+ * \brief Compute the diffusive fluxes within the different phases.
+ */
+ Scalar computeDiffusivePhaseFlux(const FluxVariables &fluxVars) const
+ {
+ // add diffusive flux of gas component in liquid phase
+ Scalar diffusivePhaseFlux = fluxVars.moleFractionGrad(phaseIdx)*fluxVars.face().normal;
+
+ return -1.0
+ * fluxVars.porousDiffCoeff(phaseIdx)
+ * fluxVars.molarDensity(phaseIdx)
+ * diffusivePhaseFlux
+ * FluidSystem::molarMass(compIdx);
+ }
+
+ /*!
+ * \brief Set the Dirichlet-like conditions for the coupling
+ * and replace the existing residual
+ *
+ * \param idx vertex index for the coupling condition
+ */
+ void evalCouplingVertex_(const int scvIdx)
+ {
+ const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
+
+ // set pressure as part of the momentum coupling
+ if (this->bcTypes_(scvIdx).isCouplingOutflow(massBalanceIdx))
+ this->residual_[scvIdx][massBalanceIdx] = volVars.pressure(nPhaseIdx);
+
+ // set mass fraction; TODO: this is fixed to contiWEqIdx so far!
+ if (this->bcTypes_(scvIdx).isCouplingOutflow(contiWEqIdx))
+ this->residual_[scvIdx][contiWEqIdx] = volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
+ }
+
+ /*!
+ * \brief Check if one of the boundary conditions is coupling.
+ */
+ bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isCouplingInflow(eqIdx) || bcTypes.isCouplingOutflow(eqIdx))
+ return true;
+ return false;
+ }
+
+ /*!
+ * \brief Check if one of the boundary conditions is mortar coupling.
+ */
+ bool boundaryHasMortarCoupling_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isMortarCoupling(eqIdx))
+ return true;
+ return false;
+ }
+
+ /*!
+ * \brief Check if one of the boundary conditions is Neumann.
+ */
+ bool boundaryHasNeumann_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isNeumann(eqIdx))
+ return true;
+ }
+ return false;
+ }
+
+ Scalar elementFluxes(const int scvIdx)
+ {
+ return elementFluxes_[scvIdx];
+ }
+
+protected:
+ ElementFluxVector elementFluxes_;
+};
+
+}
+
+#endif
diff --git a/dumux/multidomain/couplinglocalresiduals/Makefile.am b/dumux/multidomain/couplinglocalresiduals/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..9e02ea7ece61d2111c2f16c3e6c7c711946f70a3
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/Makefile.am
@@ -0,0 +1,4 @@
+couplinglocalresidualsdir = $(includedir)/dumux/multidomain/couplinglocalresiduals
+couplinglocalresiduals_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bc1fe7748b6e4fe7964d90ab8f176d764a49fceb
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
@@ -0,0 +1,243 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Klaus Mosthaf *
+ * Copyright (C) 2011 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Calculates the residual of models based on the box scheme element-wise.
+ */
+#ifndef DUMUX_BOX_COUPLING_LOCAL_RESIDUAL_HH
+#define DUMUX_BOX_COUPLING_LOCAL_RESIDUAL_HH
+
+#include <dumux/implicit/box/boxlocalresidual.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxModel
+ * \ingroup BoxLocalResidual
+ * \brief Element-wise calculation of the residual matrix for models
+ * based on the box scheme .
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class BoxCouplingLocalResidual : public BoxLocalResidual<TypeTag>
+{
+private:
+ typedef BoxCouplingLocalResidual<TypeTag> ThisType;
+ typedef BoxLocalResidual<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::Grid::ctype CoordScalar;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
+
+ typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
+ typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
+
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename Element::Geometry Geometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
+ typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+
+ // copying the local residual class is not a good idea
+ BoxCouplingLocalResidual(const BoxCouplingLocalResidual &);
+
+public:
+ BoxCouplingLocalResidual()
+ { }
+
+ ~BoxCouplingLocalResidual()
+ { }
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero. Gets a solution vector computed by PDELab
+ *
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ * \param elementSolVector The local solution for the element using PDELab ordering
+ */
+ template<typename ElemSolVectorType>
+ void evalPDELab(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElemSolVectorType& elementSolVector,
+ ElementVolumeVariables& volVarsPrev,
+ ElementVolumeVariables& volVarsCur)
+ {
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeometry;
+
+ volVarsPrev.update(this->problem_(),
+ element,
+ fvGeometry,
+ true /* oldSol? */);
+ volVarsCur.updatePDELab(this->problem_(),
+ element,
+ fvGeometry,
+ elementSolVector);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeometry);
+
+ asImp_().evalPDELab(element, fvGeometry, volVarsPrev, volVarsCur, bcTypes);
+ }
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero without taking boundary conditions into account.
+ * This is required for the flux calculation at the interface
+ * (called from the coupled problem). Calls evalPDELab with the
+ * required removal of the stabilization at the boundary (stokes).
+ *
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ */
+ void evalNoBoundary(const Element &element,
+ const FVElementGeometry fvGeometry,
+ ElementVolumeVariables& volVarsPrev,
+ ElementVolumeVariables& volVarsCur)
+ {
+ volVarsPrev.update(this->problem_(),
+ element,
+ fvGeometry,
+ true /* oldSol? */);
+ volVarsCur.update(this->problem_(),
+ element,
+ fvGeometry,
+ false /* oldSol? */);
+
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeometry);
+
+ asImp_().evalPDELab(element, fvGeometry, volVarsPrev, volVarsCur, bcTypes);
+ }
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero. Calls evalBoundaryPDELab,
+ * where the stabilization of the mass balance (stokes)
+ * is removed. No further boundary conditions are employed.
+ *
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ * \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
+ * \param curVolVars The volume averaged variables for all
+ * sub-control volumes of the element at the current
+ * time level
+ * \param bcTypes The types of the boundary conditions for all
+ * vertices of the element
+ */
+ void evalPDELab(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const ElementVolumeVariables &prevVolVars,
+ const ElementVolumeVariables &curVolVars,
+ const ElementBoundaryTypes &bcTypes)
+ {
+ const int numVerts = fvGeometry.numScv;
+#if HAVE_VALGRIND
+ for (int i=0; i < numVerts; i++) {
+ Valgrind::CheckDefined(prevVolVars[i]);
+ Valgrind::CheckDefined(curVolVars[i]);
+ }
+#endif // HAVE_VALGRIND
+
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeometry;
+ this->bcTypesPtr_ = &bcTypes;
+ this->prevVolVarsPtr_ = &prevVolVars;
+ this->curVolVarsPtr_ = &curVolVars;
+
+ // reset residual
+ this->residual_.resize(numVerts);
+ this->storageTerm_.resize(numVerts);
+
+ this->residual_ = 0;
+ this->storageTerm_ = 0;
+
+ asImp_().evalFluxes_();
+ asImp_().evalVolumeTerms_();
+
+ // evaluate the boundary (modified version)
+ asImp_().evalBoundaryPDELab_();
+
+#if HAVE_VALGRIND
+ for (int i=0; i < numVerts; i++)
+ Valgrind::CheckDefined(this->residual_[i]);
+#endif // HAVE_VALGRIND
+ }
+
+protected:
+ /*!
+ * /brief Empty method, has to be overwritten if required.
+ * Called e.g. for the removal of the stabilization of the
+ * stokes model.
+ *
+ */
+ void evalBoundaryPDELab_()
+ { }
+
+ Implementation &asImp_()
+ {
+ assert(static_cast<Implementation*>(this) != 0);
+ return *static_cast<Implementation*>(this);
+ }
+
+ const Implementation &asImp_() const
+ {
+ assert(static_cast<const Implementation*>(this) != 0);
+ return *static_cast<const Implementation*>(this);
+ }
+};
+
+}
+
+#endif