From 93ed96bfb42c0f6410c6a221c7acd59316437393 Mon Sep 17 00:00:00 2001
From: Benjamin Faigle <benjamin.faigle@posteo.de>
Date: Fri, 13 Jan 2012 14:03:02 +0000
Subject: [PATCH] corrected problem() into problem_ introduced property
 specifying if faces are regarded from both sides: compositional models do not
 yet run without

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7365 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/decoupled/2p2c/2p2cproperties.hh        |  2 +
 dumux/decoupled/common/decoupledproperties.hh |  4 ++
 dumux/decoupled/common/fv/fvpressure.hh       | 53 ++++++++++---------
 3 files changed, 35 insertions(+), 24 deletions(-)

diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index 33e186714c..6977a36ffb 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -136,6 +136,8 @@ SET_PROP(DecoupledTwoPTwoC, TransportSolutionType)
 };
 
 SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCompressibility, true);
+//! Faces are regarded from both sides
+SET_BOOL_PROP(DecoupledTwoPTwoC, VisitFacesOnlyOnce, false);
 
 SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false);
 
diff --git a/dumux/decoupled/common/decoupledproperties.hh b/dumux/decoupled/common/decoupledproperties.hh
index b7462d417e..9bd593b995 100644
--- a/dumux/decoupled/common/decoupledproperties.hh
+++ b/dumux/decoupled/common/decoupledproperties.hh
@@ -71,6 +71,7 @@ NEW_PROP_TAG( Variables); //!< The type of the container of global variables
 NEW_PROP_TAG(TimeManager);  //!< Manages the simulation time
 NEW_PROP_TAG(BoundaryTypes); //!< Stores the boundary types of a single degree of freedom
 NEW_PROP_TAG( CellData );//!< Defines data object to be stored
+NEW_PROP_TAG( VisitFacesOnlyOnce); //!< Indicates if faces are only regarded from one side
 }
 }
 
@@ -102,6 +103,9 @@ public:
 //! Disables Grid Adaptivity as standard
 SET_BOOL_PROP(DecoupledModel, AdaptiveGrid, false);
 
+//! Faces are only regarded from one side and not from both cells
+SET_BOOL_PROP(DecoupledModel, VisitFacesOnlyOnce, true);
+
 NEW_PROP_TAG(MaxIntersections);   //!< maximum number of intersections per element
 SET_PROP(DecoupledModel, MaxIntersections)
 {
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index 8352bf1f4a..303dd052fa 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -131,13 +131,13 @@ public:
     //! Function needed for restart option.
     void serializeEntity(std::ostream &outstream, const Element &element)
     {
-        int globalIdx = problem().variables().index(element);
+        int globalIdx = problem_.variables().index(element);
         outstream << pressure_[globalIdx][0];
     }
 
     void deserializeEntity(std::istream &instream, const Element &element)
     {
-        int globalIdx = problem().variables().index(element);
+        int globalIdx = problem_.variables().index(element);
         instream >> pressure_[globalIdx][0];
     }
     //@}
@@ -174,18 +174,18 @@ template<class TypeTag>
 void FVPressure<TypeTag>::initializeMatrix()
 {
     // determine matrix row sizes
-    ElementIterator eItEnd = problem().gridView().template end<0> ();
-    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+    ElementIterator eItEnd = problem_.gridView().template end<0> ();
+    for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
     {
         // cell index
-        int globalIdxI = problem().variables().index(*eIt);
+        int globalIdxI = problem_.variables().index(*eIt);
 
         // initialize row size
         int rowSize = 1;
 
         // run through all intersections with neighbors
-        IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
-        for (IntersectionIterator isIt = problem().gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+        IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+        for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
         {
             if (isIt->neighbor())
                 rowSize++;
@@ -195,22 +195,22 @@ void FVPressure<TypeTag>::initializeMatrix()
     A_.endrowsizes();
 
     // determine position of matrix entries
-    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+    for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
     {
         // cell index
-        int globalIdxI = problem().variables().index(*eIt);
+        int globalIdxI = problem_.variables().index(*eIt);
 
         // add diagonal index
         A_.addindex(globalIdxI, globalIdxI);
 
         // run through all intersections with neighbors
-        IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
-        for (IntersectionIterator isIt = problem().gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+        IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+        for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
             if (isIt->neighbor())
             {
                 // access neighbor
                 ElementPointer outside = isIt->outside();
-                int globalIdxJ = problem().variables().index(*outside);
+                int globalIdxJ = problem_.variables().index(*outside);
 
                 // add off diagonal index
                 A_.addindex(globalIdxI, globalIdxJ);
@@ -236,12 +236,12 @@ void FVPressure<TypeTag>::assemble(bool first)
     A_ = 0;
     f_ = 0;
 
-    ElementIterator eItEnd = problem().gridView().template end<0> ();
-    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+    ElementIterator eItEnd = problem_.gridView().template end<0> ();
+    for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
     {
         // cell information
-        int globalIdxI = problem().variables().index(*eIt);
-        CellData& cellDataI = problem().variables().cellData(globalIdxI);
+        int globalIdxI = problem_.variables().index(*eIt);
+        CellData& cellDataI = problem_.variables().cellData(globalIdxI);
 
         Dune::FieldVector<Scalar, 2> entries(0.);
 
@@ -251,16 +251,17 @@ void FVPressure<TypeTag>::assemble(bool first)
 
         /*****  flux term ***********/
         // iterate over all faces of the cell
-        IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
-        for (IntersectionIterator isIt = problem().gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+        IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+        for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
         {
             /************* handle interior face *****************/
             if (isIt->neighbor())
             {
-                int globalIdxJ = problem().variables().index(*(isIt->outside()));
+                int globalIdxJ = problem_.variables().index(*(isIt->outside()));
 
                 //calculate only from one side, but add matrix entries for both sides
-                if (globalIdxI > globalIdxJ)
+                if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) &&
+                        (globalIdxI > globalIdxJ))
                     continue;
 
                 entries = 0;
@@ -269,14 +270,18 @@ void FVPressure<TypeTag>::assemble(bool first)
 
                 //set right hand side
                 f_[globalIdxI] -= entries[rhs];
-                f_[globalIdxJ] += entries[rhs];
 
                 // set diagonal entry
                 A_[globalIdxI][globalIdxI] += entries[matrix];
-                A_[globalIdxJ][globalIdxJ] += entries[matrix];
                 // set off-diagonal entry
                 A_[globalIdxI][globalIdxJ] = -entries[matrix];
-                A_[globalIdxJ][globalIdxI] = -entries[matrix];
+
+                if(GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce))
+                {
+                    f_[globalIdxJ] += entries[rhs];
+                    A_[globalIdxJ][globalIdxJ] += entries[matrix];
+                    A_[globalIdxJ][globalIdxI] = -entries[matrix];
+                }
             }   // end neighbor
 
 
@@ -317,7 +322,7 @@ void FVPressure<TypeTag>::solve()
     if (verboseLevelSolver)
         std::cout << __FILE__ <<": solve for pressure" << std::endl;
 
-    Solver solver(problem());
+    Solver solver(problem_);
     solver.solve(A_, pressure_, f_);
 //                    printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
 //                    printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
-- 
GitLab