fix implicit velocity output, which counted too many surrounding cells in case...

fix implicit velocity output, which counted too many surrounding cells in case of more than one phase. This also made the function completeVelocityCalculation obsolete.

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

dumux/implicit/1p/1pmodel.hh

@@ -134,8 +134,6 @@ public:
             velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, /*phaseIdx=*/0);
         }
 
-        velocityOutput.completeVelocityCalculation(*velocity);
-
         writer.attachDofData(*p, "p", isBox);
         writer.attachDofData(*K, "K", isBox);
         if (velocityOutput.enableOutput())

dumux/implicit/1p2c/1p2cmodel.hh

@@ -175,8 +175,6 @@ public:
             velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, phaseIdx);
         }
 
-        velocityOutput.completeVelocityCalculation(*velocity);
-
         writer.attachDofData(pressure, "P", isBox);
         writer.attachDofData(delp, "delp", isBox);
         if (velocityOutput.enableOutput())

dumux/implicit/2p/2pmodel.hh

@@ -187,9 +187,6 @@ public:
             velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *elemIt, nPhaseIdx);
         }
 
-        velocityOutput.completeVelocityCalculation(*velocityW);
-        velocityOutput.completeVelocityCalculation(*velocityN);
-
         writer.attachDofData(*Sn, "Sn", isBox);
         writer.attachDofData(*Sw, "Sw", isBox);
         writer.attachDofData(*pN, "pn", isBox);

dumux/implicit/2p2c/2p2cmodel.hh

@@ -375,9 +375,6 @@ public:
 
         } // loop over elements
 
-        velocityOutput.completeVelocityCalculation(*velocityW);
-        velocityOutput.completeVelocityCalculation(*velocityN);
-
         writer.attachDofData(*sN,     "Sn", isBox);
         writer.attachDofData(*sW,     "Sw", isBox);
         writer.attachDofData(*pN,     "pN", isBox);

dumux/implicit/3p3c/3p3cmodel.hh

@@ -379,10 +379,6 @@ public:
 
         }
 
-        velocityOutput.completeVelocityCalculation(*velocityW);
-        velocityOutput.completeVelocityCalculation(*velocityN);
-        velocityOutput.completeVelocityCalculation(*velocityG);
-
         writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
         writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
         writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);

dumux/implicit/common/implicitvelocityoutput.hh

@@ -70,26 +70,41 @@ public:
      * \param problem The problem to be solved
      */
     ImplicitVelocityOutput(const Problem& problem)
-        : problem_(problem)
+    : problem_(problem)
     {
         // check, if velocity output can be used (works only for cubes so far)
         velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
+        if (velocityOutput_ && isBox)
+        {
+            cellNum_.assign(problem_.gridView().size(dofCodim), 0);
+        }
+
         ElementIterator elemIt = problem_.gridView().template begin<0>();
         ElementIterator elemEndIt = problem_.gridView().template end<0>();
         for (; elemIt != elemEndIt; ++elemIt)
         {
-            if (elemIt->geometry().type().isCube() == false){
+            if (elemIt->geometry().type().isCube() == false)
+            {
                 velocityOutput_ = false;
             }
+
+            if (isBox)
+            {
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(problem_.gridView(), *elemIt);
+
+                // transform vertex velocities from local to global coordinates
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = problem_.vertexMapper().map(*elemIt, scvIdx, dofCodim);
+
+                    cellNum_[globalIdx] += 1;
+                }
+            }
         }
 
         if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
             std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
-
-        if (velocityOutput_ && isBox)
-        {
-            cellNum_.resize(problem_.gridView().size(dofCodim), 0);
-        }
     }
 
     bool enableOutput()
@@ -137,7 +152,7 @@ public:
                                            fvGeometry,
                                            faceIdx,
                                            elemVolVars);
-                    
+
                     const GlobalPosition globalNormal = fluxVars.face().normal;
 
                     GlobalPosition localNormal(0);
@@ -149,8 +164,7 @@ public:
 
                     // Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
                     // face in the reference element.
-                    Scalar flux;
-                    flux = fluxVars.volumeFlux(phaseIdx) / localArea;
+                    Scalar flux = fluxVars.volumeFlux(phaseIdx) / localArea;
 
                     // transform the normal Darcy velocity into a vector
                     tmpVelocity = localNormal;
@@ -168,11 +182,9 @@ public:
                     Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
 
                     jacobianT2.mtv(scvVelocities[scvIdx], scvVelocity);
-                    scvVelocity /= element.geometry().integrationElement(localPos);
+                    scvVelocity /= element.geometry().integrationElement(localPos)*cellNum_[globalIdx];
                     // add up the wetting phase subcontrolvolume velocities for each vertex
                     velocity[globalIdx] += scvVelocity;
-
-                    cellNum_[globalIdx] +=1;
                 }
             }
             else
@@ -206,12 +218,12 @@ public:
                                                fvGeometry,
                                                faceIdx,
                                                elemVolVars,true);
-                        
+
                         Scalar flux = fluxVars.volumeFlux(phaseIdx);
                         scvVelocities[faceIdx] = flux;
                     }
                 }
-                
+
                 Dune::FieldVector < CoordScalar, dim > refVelocity(0);
                 for (int i = 0; i < dim; i++)
                     refVelocity[i] = 0.5 * (scvVelocities[2*i + 1] - scvVelocities[2*i]);
@@ -228,32 +240,10 @@ public:
         } // velocity output
     }
 
-    template<class VelocityVector>
-    void completeVelocityCalculation(VelocityVector& velocity)
-    {
-        if (velocityOutput_ && isBox)
-        {
-            unsigned int size =  velocity.size();
-
-            assert(size != cellNum_.size());
-
-
-            // divide the vertex velocities by the number of adjacent scvs i.e. cells
-            for(unsigned int globalIdx = 0; globalIdx < size; ++globalIdx)
-            {
-                velocity[globalIdx] /= cellNum_[globalIdx];
-            }
-            fluxVariables_.clear();
-        }
-    }
 protected:
     const Problem& problem_;
-
-    // parameters given in constructor
     bool velocityOutput_;
     std::vector<int> cellNum_;
-
-    std::vector<FluxVariables> fluxVariables_;
 };
 
 }

dumux/implicit/mpnc/mpncvtkwritercommon.hh

@@ -65,7 +65,7 @@ class MPNCVtkWriterCommon : public MPNCVtkWriterModule<TypeTag>
     typedef Dune::array<DimField, numPhases> PhaseDimField;
     enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
     enum { dofCodim = isBox ? dim : 0 };
-    
+
 public:
     MPNCVtkWriterCommon(const Problem &problem)
     : ParentType(problem), velocityOutput_(problem)
@@ -99,7 +99,7 @@ public:
 
         if (velocityOutput_.enableOutput()) {
             Scalar nDofs = this->problem_.model().numDofs();
-            for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                 velocity_[phaseIdx].resize(nDofs);
                 velocity_[phaseIdx] = 0;
             }
@@ -127,7 +127,7 @@ public:
         for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
         {
             int globalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
-            
+
             const VolumeVariables &volVars = elemVolVars[scvIdx];
 
             if (porosityOutput_) porosity_[globalIdx] = volVars.porosity();
@@ -203,7 +203,6 @@ public:
 
         if (velocityOutput_.enableOutput()) {
             for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                velocityOutput_.completeVelocityCalculation(velocity_[phaseIdx]);
                 std::ostringstream oss;
                 oss << "velocity_" << FluidSystem::phaseName(phaseIdx);
                 writer.attachDofData(velocity_[phaseIdx],

dumux/implicit/richards/richardsmodel.hh

@@ -206,8 +206,6 @@ public:
             velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, /*phaseIdx=*/0);
         }
 
-        velocityOutput.completeVelocityCalculation(*velocity);
-
         writer.attachDofData(*Sn, "Sn", isBox);
         writer.attachDofData(*Sw, "Sw", isBox);
         writer.attachDofData(*pN, "pn", isBox);