diff --git a/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..229205486e6fc5b8e28c99805796307d8a83b893
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
@@ -0,0 +1,344 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Klaus Mosthaf *
+ * Copyright (C) 2008 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 Supplements the TwoPTwoCNILocalResidual by the required functions for a coupled application.
+ *
+ */
+#ifndef DUMUX_2P2CNI_COUPLING_LOCAL_RESIDUAL_HH
+#define DUMUX_2P2CNI_COUPLING_LOCAL_RESIDUAL_HH
+
+#include <dumux/implicit/2p2cni/2p2cnilocalresidual.hh>
+
+#include <dumux/implicit/2p2cnicoupling/2p2cnicouplingproperties.hh>
+
+#define VELOCITY_OUTPUT 1 // uncomment this line if an output of the velocity is needed
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief Supplements the TwoPTwoCNIModel by the required functions for a coupled application.
+ */
+template<class TypeTag>
+class TwoPTwoCNICouplingLocalResidual : public TwoPTwoCNILocalResidual<TypeTag>
+{
+ typedef TwoPTwoCNICouplingLocalResidual<TypeTag> ThisType;
+ typedef TwoPTwoCNILocalResidual<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,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+ enum {
+ massBalanceIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
+ contiWEqIdx = Indices::contiWEqIdx,
+ energyEqIdx = Indices::energyEqIdx
+ };
+ 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.
+ */
+ TwoPTwoCNICouplingLocalResidual()
+ { };
+
+// void evalBoundary_()
+// {
+// ParentType::evalBoundary_();
+//
+// // evaluate Dirichlet-like coupling conditions
+// for (int idx = 0; idx < this->fvGeometry_().numVertices; idx++)
+// if (boundaryHasCoupling_(this->bcTypes_(idx)))
+// evalCouplingVertex_(idx);
+// }
+
+ 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 elemVertIdx = refElement.subEntity(faceIdx, 1, faceVertIdx, dim);
+ const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ // 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->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 a Neumann condition is set
+ //in case of Neumann Null the boundary flux must not be calculated
+ if(this->model_().onBoundary(this->element_(), i) && i!=elemVertIdx)
+ if (!this->bcTypes_(i).isNeumann(equationIdx) && !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;
+ }
+
+
+ Scalar computePhaseStorage(const int scvIdx, bool usePrevSol) const
+ {
+ const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+ // compute storage term of all components within all phases
+ Scalar phaseStorage = volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFraction(phaseIdx, compIdx);
+ phaseStorage *= volVars.porosity();
+
+ return phaseStorage;
+ }
+
+ 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;
+ }
+ }
+
+ 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;
+ }
+
+ Scalar computeDiffusivePhaseFlux(const FluxVariables &fluxVars) const
+ {
+ // add diffusive flux of gas component in liquid phase
+ Scalar diffusivePhaseFlux = fluxVars.moleFractionGrad(phaseIdx)*fluxVars.face().normal;
+ diffusivePhaseFlux *= -fluxVars.porousDiffCoeff(phaseIdx) *
+ fluxVars.molarDensity(phaseIdx);
+ diffusivePhaseFlux *= FluidSystem::molarMass(compIdx);
+
+ return diffusivePhaseFlux;
+ }
+
+ /*!
+ * \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);
+
+ // set temperature as part of the coupling
+ if (this->bcTypes_(scvIdx).isCouplingOutflow(energyEqIdx))
+ this->residual_[scvIdx][energyEqIdx] = volVars.temperature();
+ }
+
+ /*!
+ * \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 Neumann.
+ */
+ bool boundaryHasNeumann_(const BoundaryTypes& bcTypes) const
+ {
+ bool hasNeumann = false;
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isNeumann(eqIdx) || bcTypes.isNeumann(eqIdx))
+ hasNeumann = true;
+ return hasNeumann;
+ }
+
+ Scalar elementFluxes(const int scvIdx)
+ {
+ return elementFluxes_[scvIdx];
+ }
+
+private:
+ ElementFluxVector elementFluxes_;
+};
+
+}
+
+#endif
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..db4f3ab25c23e73777273bde4f2d7efe8f69b494
--- /dev/null
+++ b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
@@ -0,0 +1,418 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2012 by Katherina Baber, Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
+ * 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the stokes box model.
+ */
+
+#ifndef DUMUX_STOKESNCNI_COUPLING_LOCAL_RESIDUAL_HH
+#define DUMUX_STOKESNCNI_COUPLING_LOCAL_RESIDUAL_HH
+
+#include <dumux/freeflow/stokesncni/stokesncnilocalresidual.hh>
+#include <dumux/freeflow/stokesncni/stokesncnimodel.hh>
+
+namespace Dumux
+{
+ /*!
+ * \ingroup BoxStokesModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the Stokes box model.
+ *
+ * This class is also used for the stokes transport
+ * model, which means that it uses static polymorphism.
+ */
+ template<class TypeTag>
+ class StokesncniCouplingLocalResidual : public StokesncniLocalResidual<TypeTag>
+ {
+ protected:
+ typedef StokesncniCouplingLocalResidual<TypeTag> ThisType;
+ typedef StokesncniLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ numComponents = FluidSystem::numComponents
+ };
+ enum {
+ //indices of the equations
+ massBalanceIdx = Indices::massBalanceIdx, //!< Index of the mass balance
+
+ momentumXIdx = Indices::momentumXIdx, //!< Index of the x-component of the momentum balance
+ momentumYIdx = Indices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ momentumZIdx = Indices::momentumZIdx, //!< Index of the z-component of the momentum balance
+ lastMomentumIdx = Indices::lastMomentumIdx, //!< Index of the last component of the momentum balance
+ transportEqIdx = Indices::transportEqIdx,//!< Index of the transport equation
+ energyEqIdx = Indices::energyEqIdx
+ };
+ enum {
+ //indices of phase and transported component
+ phaseIdx = Indices::phaseIdx,
+ transportCompIdx = Indices::transportCompIdx,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum {
+ dimXIdx = Indices::dimXIdx, //!< Index for the first component of a vector
+ dimYIdx = Indices::dimYIdx, //!< Index for the second component of a vector
+ dimZIdx = Indices::dimZIdx //!< Index for the third component of a vector
+ };
+
+ 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 GridView::ctype CoordScalar;
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GridView::Intersection Intersection;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+ public:
+ // the stokes model needs a modified treatment of the BCs
+ void evalBoundary_()
+ {
+ assert(this->residual_.size() == this->fvGeometry_().numScv);
+ const ReferenceElement &refElement = ReferenceElements::general(this->element_().geometry().type());
+
+ // loop over vertices of the element
+ for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
+ {
+ // consider only SCVs on the boundary
+ if (this->fvGeometry_().subContVol[idx].inner)
+ continue;
+
+ // important at corners of the grid
+ DimVector momentumResidual(0.0);
+ DimVector averagedNormal(0.0);
+ int numberOfOuterFaces = 0;
+ // evaluate boundary conditions for the intersections of
+ // the current element
+ const BoundaryTypes &bcTypes = this->bcTypes_(idx);
+ 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)
+ {
+ // only evaluate, if we consider the same face vertex as in the outer
+ // loop over the element vertices
+ if (refElement.subEntity(faceIdx, 1, faceVertIdx, dim)
+ != idx)
+ continue;
+
+ const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ const FluxVariables boundaryVars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ boundaryFaceIdx,
+ this->curVolVars_(),
+ true);
+
+ // the computed residual of the momentum equations is stored
+ // into momentumResidual for the replacement of the mass balance
+ // in case of Dirichlet conditions for the momentum balance;
+ // the fluxes at the boundary are added in the second step
+ if (this->momentumBalanceDirichlet_(bcTypes))
+ {
+ DimVector muGradVelNormal(0.);
+ const DimVector &boundaryFaceNormal =
+ boundaryVars.face().normal;
+
+ boundaryVars.velocityGrad().umv(boundaryFaceNormal, muGradVelNormal);
+ muGradVelNormal *= boundaryVars.dynamicViscosity();
+
+ for (int i=0; i < this->residual_.size(); i++)
+ Valgrind::CheckDefined(this->residual_[i]);
+ for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ momentumResidual[dimIdx] = this->residual_[idx][momentumXIdx+dimIdx];
+
+ //Sign is right!!!: boundary flux: -mu grad v n
+ //but to compensate outernormal -> residual - (-mu grad v n)
+ momentumResidual += muGradVelNormal;
+ averagedNormal += boundaryFaceNormal;
+ }
+
+ // evaluate fluxes at a single boundary segment
+ asImp_()->evalNeumannSegment_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ asImp_()->evalOutflowSegment_(isIt, idx, boundaryFaceIdx, boundaryVars);
+
+ //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->momentumBalanceHasNeumann_(this->bcTypes_(idx)))
+ {
+ const GlobalPosition& globalPos = this->fvGeometry_().subContVol[idx].global;
+ //problem specific function, in problem orientation of interface is known
+ if(this->problem_().isInterfaceCornerPoint(globalPos))
+ {
+ PrimaryVariables priVars(0.0);
+ // DimVector faceCoord = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
+ // std::cout<<faceCoord<<std::endl;
+ //calculate the actual boundary fluxes and add to residual (only for momentum equation, mass balance already has outflow)
+ asImp_()->computeFlux(priVars, boundaryFaceIdx, true/*on boundary*/);
+ for(int equationIdx = 0; equationIdx < numEq; ++equationIdx)
+ {
+ if(equationIdx == massBalanceIdx)
+ continue;
+ this->residual_[idx][equationIdx] += priVars[equationIdx];
+ }
+ }
+ }
+ // Beavers-Joseph condition at the coupling boundary/interface
+ if(boundaryHasCoupling_(bcTypes))
+ {
+ evalBeaversJoseph_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ asImp_()->evalCouplingVertex_(isIt, idx, boundaryFaceIdx, boundaryVars);
+ }
+
+ // count the number of outer faces to determine, if we are on
+ // a corner point and if an interpolation should be done
+ numberOfOuterFaces++;
+ } // end loop over face vertices
+ } // end loop over intersections
+
+ // replace defect at the corner points of the grid
+ // by the interpolation of the primary variables
+ if(!bcTypes.isDirichlet(massBalanceIdx))
+ {
+ if (this->momentumBalanceDirichlet_(bcTypes))
+ this->replaceMassbalanceResidual_(momentumResidual, averagedNormal, idx);
+ else // de-stabilize (remove alpha*grad p - alpha div f
+ // from computeFlux on the boundary)
+ this->removeStabilizationAtBoundary_(idx);
+ }
+ if (numberOfOuterFaces == 2)
+ this->interpolateCornerPoints_(bcTypes, idx);
+ } // end loop over element vertices
+
+ // evaluate the dirichlet conditions of the element
+ if (this->bcTypes_().hasDirichlet())
+ asImp_()->evalDirichlet_();
+ }
+
+ void evalBoundaryPDELab_()
+ {
+ // loop over vertices of the element
+ for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
+ {
+ // consider only SCVs on the boundary
+ if (this->fvGeometry_().subContVol[idx].inner)
+ continue;
+
+ this->removeStabilizationAtBoundary_(idx);
+ } // end loop over vertices
+ }
+
+ protected:
+ /*!
+ * \brief Evaluate one part of the Dirichlet-like coupling conditions for a single
+ * sub-control volume face; rest is done in the local coupling operator
+ */
+ void evalCouplingVertex_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars)
+ {
+ const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ // TODO: beaversJosephCoeff is used to determine, if we are on a coupling face
+ // this is important at the corners. However, a better way should be found.
+
+ // check if BJ-coeff is not zero to decide, if coupling condition
+ // for the momentum balance (Dirichlet vor v.n) has to be set;
+ // may be important at corners
+ const Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+ // set velocity normal to the interface
+ if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
+ this->residual_[scvIdx][momentumYIdx] =
+ volVars.velocity() *
+ boundaryVars.face().normal /
+ boundaryVars.face().normal.two_norm();
+ Valgrind::CheckDefined(this->residual_[scvIdx][momentumYIdx]);
+
+ // add pressure correction - required for pressure coupling,
+ // if p.n comes from the pm
+ if (bcTypes.isCouplingOutflow(momentumYIdx) && beaversJosephCoeff)
+ {
+ DimVector pressureCorrection(isIt->centerUnitOuterNormal());
+ pressureCorrection *= volVars.pressure(); // TODO: 3D
+ this->residual_[scvIdx][momentumYIdx] += pressureCorrection[momentumYIdx]*
+ boundaryVars.face().area;
+ }
+
+ // set mole fraction for the transported components
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ int eqIdx = dim + compIdx;
+ if (eqIdx != massBalanceIdx) {
+ if (bcTypes.isCouplingOutflow(eqIdx))
+ {
+ if(useMoles)
+ this->residual_[scvIdx][eqIdx] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
+ else
+ this->residual_[scvIdx][eqIdx] = volVars.fluidState().massFraction(phaseIdx, compIdx);
+ Valgrind::CheckDefined(this->residual_[scvIdx][eqIdx]);
+ }
+ }
+ }
+ // set temperature
+ if (bcTypes.isCouplingOutflow(energyEqIdx))
+ {
+ this->residual_[scvIdx][energyEqIdx] = volVars.temperature();
+ Valgrind::CheckDefined(this->residual_[scvIdx][energyEqIdx]);
+ }
+ }
+
+ void evalBeaversJoseph_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars) //TODO: required
+ {
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeff(this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+
+ // only enter here, if we are on a coupling boundary (see top)
+ // and the BJ coefficient is not zero
+ if (beaversJosephCoeff)
+ {
+ const Scalar Kxx = this->problem_().permeability(this->element_(),
+ this->fvGeometry_(),
+ scvIdx);
+ beaversJosephCoeff /= sqrt(Kxx);
+ const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
+
+ // implementation as NEUMANN condition /////////////////////////////////////////////
+ // (v.n)n
+ if (bcTypes.isCouplingOutflow(momentumXIdx))
+ {
+ const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
+ DimVector normalV = elementUnitNormal;
+ normalV *= normalComp; // v*n*n
+
+ // v_tau = v - (v.n)n
+ const DimVector tangentialV = boundaryVars.velocity() - normalV;
+ const Scalar boundaryFaceArea = boundaryVars.face().area;
+
+ for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
+ boundaryFaceArea *
+ tangentialV[dimIdx] *
+ boundaryVars.dynamicViscosity();
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //normal component has only to be set if no coupling conditions are defined
+ //set NEUMANN flux (set equal to pressure in problem)
+ // PrimaryVariables priVars(0.0);
+ // this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
+ // *isIt, scvIdx, boundaryFaceIdx);
+ // for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ // this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
+ // boundaryFaceArea;
+ }
+ if (bcTypes.isCouplingInflow(momentumXIdx)) //TODO: 3D
+ {
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ //IMPLEMENTATION AS DIRICHLET CONDITION
+ //tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
+ DimVector tangentialVelGrad(0);
+ boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
+ tangentialVelGrad /= -beaversJosephCoeff; // was - before
+ this->residual_[scvIdx][momentumXIdx] =
+ tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //for testing Beavers and Joseph without adjacent porous medium set vy to 0
+ // Scalar normalVel(0);
+ // this->residual_[scvIdx][momentumYIdx] = this->curPrimaryVars_(scvIdx)[momentumYIdx] - normalVel;
+ ////////////////////////////////////////////////////////////////////////////////////
+ }
+ }
+ }
+
+ // return true, if at least one equation on the boundary has a coupling condition
+ bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isCouplingInflow(eqIdx) || bcTypes.isCouplingOutflow(eqIdx))
+ return true;
+ return false;
+ }
+
+ private:
+ Implementation *asImp_()
+ { return static_cast<Implementation *>(this); }
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *>(this); }
+
+ };
+
+}
+
+#endif