implicit branch: fix outflow conditions for cell-centered disc. Unified 1p2c test works now.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9774 2fb0f335-1f38-0410-981e-8018bf24f1b0

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;
             }
         }
     }

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);
             

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;
     }
 
     /*!