From fae3859c7be9ce9547db734e5488e79c23269cfa Mon Sep 17 00:00:00 2001 From: Bernd Flemisch <bernd@iws.uni-stuttgart.de> Date: Wed, 22 May 2013 11:31:20 +0000 Subject: [PATCH] implicit mpnc: set cell-centered Dirichlet conditions in the strong sense, until someone manages to do them in a correct way git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10729 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- dumux/implicit/mpnc/mpnclocalresidual.hh | 57 ++++++++++++++++++++---- 1 file changed, 49 insertions(+), 8 deletions(-) diff --git a/dumux/implicit/mpnc/mpnclocalresidual.hh b/dumux/implicit/mpnc/mpnclocalresidual.hh index 3ce957fd1f..b66685bba9 100644 --- a/dumux/implicit/mpnc/mpnclocalresidual.hh +++ b/dumux/implicit/mpnc/mpnclocalresidual.hh @@ -53,18 +53,21 @@ protected: typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)}; + enum {numEq = GET_PROP_VALUE(TypeTag, NumEq)}; enum {enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy)}; enum {enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy)}; enum {enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic)}; enum {phase0NcpIdx = Indices::phase0NcpIdx}; typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::IntersectionIterator IntersectionIterator; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; 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 GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes; + typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef MPNCLocalResidualEnergy<TypeTag, enableEnergy, enableKineticEnergy> EnergyResid; typedef MPNCLocalResidualMass<TypeTag, enableKinetic> MassResid; @@ -225,20 +228,58 @@ public: curVolVars, bcType); - for (int i = 0; i < this->fvGeometry_().numScv; ++i) { - // add the two auxiliary equations, make sure that the - // dirichlet boundary condition is conserved - for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - { - if (!bcType[i].isDirichlet(phase0NcpIdx + phaseIdx)) + if (GET_PROP_VALUE(TypeTag, ImplicitIsBox) + || !bcType.hasDirichlet()) + { + for (int i = 0; i < this->fvGeometry_().numScv; ++i) { + // add the two auxiliary equations, make sure that the + // dirichlet boundary condition is conserved + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - this->residual_[i][phase0NcpIdx + phaseIdx] = - this->curVolVars_(i).phaseNcp(phaseIdx); + if (!bcType[i].isDirichlet(phase0NcpIdx + phaseIdx)) + { + this->residual_[i][phase0NcpIdx + phaseIdx] = + this->curVolVars_(i).phaseNcp(phaseIdx); + } } } } } + + /*! + * \brief Add Dirichlet boundary conditions for a single intersection + * + * Sets the Dirichlet conditions in a strong sense, in contrast to + * the general handling in CCLocalResidual. + */ + void evalDirichletSegment_(const IntersectionIterator &isIt, + const BoundaryTypes &bcTypes) + { + // temporary vector to store the Dirichlet boundary fluxes + PrimaryVariables values; + Valgrind::SetUndefined(values); + this->problem_().dirichlet(values, *isIt); + Valgrind::CheckDefined(values); + + // set Dirichlet conditions in a strong sense + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isDirichlet(eqIdx)) + { + int pvIdx = bcTypes.eqToDirichletIndex(eqIdx); + this->residual_[0][eqIdx] + = this->curPriVar_(0, pvIdx) - values[pvIdx]; + } + else if (eqIdx >= phase0NcpIdx) + { + int phaseIdx = eqIdx - phase0NcpIdx; + this->residual_[0][eqIdx] = + this->curVolVars_(0).phaseNcp(phaseIdx); + } + } + } + protected: Implementation &asImp_() { return *static_cast<Implementation *>(this); } -- GitLab