diff --git a/dumux/implicit/cellcentered/ccfvelementgeometry.hh b/dumux/implicit/cellcentered/ccfvelementgeometry.hh index 79e88ddee03724039274a614596d16d41ff4c3f1..719f0c675ff06d4165005cc331fc63de489b7ab5 100644 --- a/dumux/implicit/cellcentered/ccfvelementgeometry.hh +++ b/dumux/implicit/cellcentered/ccfvelementgeometry.hh @@ -96,6 +96,7 @@ public: Dune::FieldVector<Vector, maxNFAP> grad; //!< derivatives of shape functions at ip Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points + unsigned faceIdx; }; typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef @@ -183,6 +184,8 @@ public: subContVolFace[k].fapIndices[0] = subContVolFace[k].i; subContVolFace[k].fapIndices[1] = subContVolFace[k].j; + + subContVolFace[k].faceIdx = it->indexInInside(); } // boundary cvf data @@ -199,7 +202,7 @@ public: GlobalPosition distVec = elementGlobal; distVec -= boundaryFace[bfIdx].ipGlobal; - distVec /= distVec.two_norm2(); + distVec /= 2.0*distVec.two_norm2(); // gradients using a two-point flux approximation for (int idx = 0; idx < 2; idx++) @@ -207,7 +210,41 @@ public: boundaryFace[bfIdx].grad[idx] = distVec; boundaryFace[bfIdx].shapeValue[idx] = 0.5; } - boundaryFace[bfIdx].grad[1] *= -1.0; + boundaryFace[bfIdx].grad[0] *= -1.0; + + boundaryFace[bfIdx].fapIndices[0] = boundaryFace[bfIdx].i; + boundaryFace[bfIdx].fapIndices[1] = boundaryFace[bfIdx].j; + } + } + + for (int nIdx = 0; nIdx < numNeighbors-1; nIdx++) + { + switch (subContVolFace[nIdx].faceIdx) + { + case 0: + boundaryFace[1].j = nIdx+1; + boundaryFace[1].fapIndices[1] = boundaryFace[1].j; + break; + case 1: + boundaryFace[0].j = nIdx+1; + boundaryFace[0].fapIndices[1] = boundaryFace[0].j; + break; + case 2: + boundaryFace[3].j = nIdx+1; + boundaryFace[3].fapIndices[1] = boundaryFace[3].j; + break; + case 3: + boundaryFace[2].j = nIdx+1; + boundaryFace[2].fapIndices[1] = boundaryFace[2].j; + break; + case 4: + boundaryFace[5].j = nIdx+1; + boundaryFace[5].fapIndices[1] = boundaryFace[5].j; + break; + case 5: + boundaryFace[4].j = nIdx+1; + boundaryFace[4].fapIndices[1] = boundaryFace[4].j; + break; } } } diff --git a/dumux/implicit/cellcentered/cclocalresidual.hh b/dumux/implicit/cellcentered/cclocalresidual.hh index 0a89113a87a88568e0fdafb9c543096cc6f20464..a7eb533be15f338e5949b586d338e0e25ce96d66 100644 --- a/dumux/implicit/cellcentered/cclocalresidual.hh +++ b/dumux/implicit/cellcentered/cclocalresidual.hh @@ -180,9 +180,10 @@ protected: // deal with outflow boundaries if (bcTypes.hasOutflow()) { + unsigned bfIdx = isIt->indexInInside(); PrimaryVariables values(0.0); this->asImp_().computeFlux(values, - /*boundaryFaceIdx=*/isIt->indexInInside(), + bfIdx, true); Valgrind::CheckDefined(values); diff --git a/test/implicit/1p2c/1p2coutflowproblem.hh b/test/implicit/1p2c/1p2coutflowproblem.hh index 195b6ebb5bfe953fe55a221b26cbec8b636c0c94..ea0239884b87a64c5ef535b345f07e3940dfbcec 100644 --- a/test/implicit/1p2c/1p2coutflowproblem.hh +++ b/test/implicit/1p2c/1p2coutflowproblem.hh @@ -51,7 +51,7 @@ NEW_TYPE_TAG(OnePTwoCOutflowCCProblem, INHERITS_FROM(CCModel, OnePTwoCOutflowPro // Set the grid type SET_PROP(OnePTwoCOutflowProblem, Grid) { -#if 0 //HAVE_UG +#if HAVE_UG typedef Dune::UGGrid<2> type; #else typedef Dune::SGrid<2, 2> type; @@ -205,17 +205,15 @@ public: if(globalPos[0] < eps_ || globalPos[0] > this->bboxMax()[0] - eps_) { values.setAllDirichlet(); - std::cout << "set Dirichlet values at pos " << globalPos << std::endl; } else { values.setAllNeumann(); - std::cout << "set Neumann values at pos " << globalPos << std::endl; } - //outflow condition for the transport equation at right boundary - //if(globalPos[0] > this->bboxMax()[0] - eps_) - // values.setOutflow(transportEqIdx); + // outflow condition for the transport equation at right boundary + if(globalPos[0] > this->bboxMax()[0] - eps_) + values.setOutflow(transportEqIdx); } /*! @@ -229,12 +227,11 @@ public: */ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { - std::cout << "set Dirichlet values " << values << " at pos " << globalPos << std::endl; initial_(values, globalPos); - std::cout << "set Dirichlet values " << values << " at pos " << globalPos << std::endl; + //condition for the N2 molefraction at left boundary - //if (globalPos[0] < 0.5*(this->bboxMin()[0] + this->bboxMax()[0])) - // values[massOrMoleFracIdx] = 2.0e-5; + if (globalPos[0] < 0.5*(this->bboxMin()[0] + this->bboxMax()[0])) + values[massOrMoleFracIdx] = 2.0e-5; } /*!