diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4e62244593cfa598f42f40a104f595bd6122c2ed
--- /dev/null
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -0,0 +1,331 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2010 by Benjamin Faigle *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_FVPRESSURE_HH
+#define DUMUX_FVPRESSURE_HH
+
+// dune environent:
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+// dumux environment
+#include "dumux/common/math.hh"
+#include <dumux/decoupled/common/impetproperties.hh>
+
+/**
+ * @file
+ * @brief Finite Volume Diffusion Model
+ * @author Benjamin Faigle, Bernd Flemisch, Jochen Fritz, Markus Wolff
+ */
+
+namespace Dumux
+{
+//! The finite volume base class for the solution of a pressure equation
+/*! \ingroup multiphase
+ * TODO:
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag> class FVPressure
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes))::ScalarSolution ScalarSolution;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
+ contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
+ };
+
+ // typedefs to abbreviate several dune classes...
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Intersection Intersection;
+
+ // convenience shortcuts for Vectors/Matrices
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
+ typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
+ // the typenames used for the stiffness matrix and solution vector
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) RHSVector;
+
+protected:
+ //initializes the matrix to store the system of equations
+ void initializeMatrix();
+
+ //solves the system of equations to get the spatial distribution of the pressure
+ void solve();
+
+ /*! \name Access functions for protected variables */
+ //@{
+ Problem& problem()
+ {
+ return problem_;
+ }
+ const Problem& problem() const
+ {
+ return problem_;
+ }
+
+ ScalarSolution& pressure()
+ { return pressure_; }
+ const ScalarSolution& pressure() const
+ { return pressure_; }
+ //@}
+public:
+ //! Public acess function for the primary variable pressure
+ const Scalar pressure(int globalIdx) const
+ { return pressure_[globalIdx]; }
+
+ //function which assembles the system of equations to be solved
+ void assemble(bool first);
+
+ void getSource(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
+
+ void getStorage(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
+
+ void getFlux(Dune::FieldVector<Scalar, 2>&, const Intersection&, const CellData&, const bool);
+
+ void getFluxOnBoundary(Dune::FieldVector<Scalar, 2>&,
+ const Intersection&, const CellData&, const bool);
+
+ //pressure solution routine: update estimate for secants, assemble, solve.
+ void update(bool solveTwice = true)
+ {
+ problem().pressureModel().assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ problem().pressureModel().solve();
+
+ return;
+ }
+
+ /*! \name general methods for serialization, output */
+ //@{
+ // serialization methods
+ //! Function needed for restart option.
+ void serializeEntity(std::ostream &outstream, const Element &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);
+ instream >> pressure_[globalIdx][0];
+ }
+ //@}
+
+ //! Constructs a FVPressure object
+ /**
+ * \param problem a problem class object
+ */
+ FVPressure(Problem& problem) :
+ problem_(problem), A_(problem.gridView().size(0), problem.gridView().size(0), (2 * dim + 1)
+ * problem.gridView().size(0), Matrix::random), f_(problem.gridView().size(0))
+ {
+ pressure_.resize(problem.gridView().size(0));
+ }
+
+protected:
+ Problem& problem_;
+ Matrix A_;
+ RHSVector f_;
+ ScalarSolution pressure_;
+
+ std::string solverName_;
+ std::string preconditionerName_;
+
+ static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+};
+
+//! initializes the matrix to store the system of equations
+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)
+ {
+ // cell index
+ 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)
+ {
+ if (isIt->neighbor())
+ rowSize++;
+ }
+ A_.setrowsize(globalIdxI, rowSize);
+ }
+ A_.endrowsizes();
+
+ // determine position of matrix entries
+ for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ 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)
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer outside = isIt->outside();
+ int globalIdxJ = problem().variables().index(*outside);
+
+ // add off diagonal index
+ A_.addindex(globalIdxI, globalIdxJ);
+ }
+ }
+ A_.endindices();
+
+ return;
+}
+
+//! function which assembles the system of equations to be solved
+/* This function assembles the Matrix and the RHS vectors to solve for
+ * a pressure field with an Finite-Volume Discretization in an implicit
+ * fasion. In the implementations to this base class, the methods
+ * getSource(), getStorage(), getFlux() and getFluxOnBoundary() have to
+ * be provided.
+ * \param first Flag if pressure field is unknown
+ */
+template<class TypeTag>
+void FVPressure<TypeTag>::assemble(bool first)
+{
+ // initialization: set matrix A_ to zero
+ A_ = 0;
+ f_ = 0;
+
+ 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);
+
+ Dune::FieldVector<Scalar, 2> entries(0.);
+
+ /***** source term ***********/
+ problem().pressureModel().getSource(entries,*eIt, cellDataI, first);
+ f_[globalIdxI] = entries[1];
+
+ /***** 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)
+ {
+ /************* handle interior face *****************/
+ if (isIt->neighbor())
+ {
+ int globalIdxJ = problem().variables().index(*(isIt->outside()));
+ problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
+
+
+ //set right hand side
+ f_[globalIdxI] -= entries[1];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[0];
+ // set off-diagonal entry
+ A_[globalIdxI][globalIdxJ] = -entries[0];
+ } // end neighbor
+
+
+ /************* boundary face ************************/
+ else
+ {
+ problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
+
+ //set right hand side
+ f_[globalIdxI] += entries[1];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[0];
+ }
+ } //end interfaces loop
+// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
+
+ /***** storage term ***********/
+ problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
+ f_[globalIdxI] += entries[1];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[0];
+ } // end grid traversal
+// printmatrix(std::cout, A_, "global stiffness matrix after assempling", "row", 11,3);
+// printvector(std::cout, f_, "right hand side", "row", 10);
+ return;
+}
+
+//! solves the system of equations to get the spatial distribution of the pressure
+template<class TypeTag>
+void FVPressure<TypeTag>::solve()
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LinearSolver)) Solver;
+
+ int verboseLevelSolver = GET_PARAM(TypeTag, int, LinearSolver, Verbosity);
+
+ if (verboseLevelSolver)
+ std::cout << __FILE__ <<": solve for pressure" << std::endl;
+
+ 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);
+// printvector(std::cout, pressure_, "pressure", "row", 200, 1, 3);
+
+ return;
+}
+
+}//end namespace Dumux
+#endif
diff --git a/dumux/decoupled/common/gridadapt.hh b/dumux/decoupled/common/gridadapt.hh
index ba99fe2f0147fd066beb415a8babe3871c9dc795..60f0c57fdd019c6914c0ea699d1eaa4a0e202738 100644
--- a/dumux/decoupled/common/gridadapt.hh
+++ b/dumux/decoupled/common/gridadapt.hh
@@ -155,9 +155,9 @@ public:
* @param indicator Vector where the refinement indicator is stored
* @param refineThreshold lower threshold where to refine
* @param coarsenThreshold upper threshold where to coarsen
- * @return
+ * @return Total ammount of marked cells
*/
- int markElements(ScalarSolutionType &indicator, const double refineThreshold, const double coarsenThreshold)
+ int markElements(const ScalarSolutionType &indicator, const double refineThreshold, const double coarsenThreshold)
{
for (LeafIterator eIt = problem_.gridView().template begin<0>();
eIt!=problem_.gridView().template end<0>(); ++eIt)
diff --git a/dumux/decoupled/common/impet.hh b/dumux/decoupled/common/impet.hh
index 42a8895bfdc1ab65419f72c02130a1da5b9e3188..6e843db2add95d7939dd4ba2681a91643db15acb 100644
--- a/dumux/decoupled/common/impet.hh
+++ b/dumux/decoupled/common/impet.hh
@@ -67,8 +67,7 @@ template<class TypeTag> class IMPET
enum
{
- dim = GridView::dimension, dimWorld = GridView::dimensionworld,
- transportEqIdx = Indices::transportEqIdx
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
@@ -99,42 +98,44 @@ public:
* \param dt time step size
* \param updateVec vector for the update values
*
- * Calculates the new pressure and velocity and determines the time step size and the update of the transported quantity for the explicit time step
+ * Calculates the new pressure and velocity and determines the time step size
+ * and the update of the transported quantity for the explicit time step.
*/
void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
{
- // the method is valid for any transported quantity.
- int transSize = problem.variables().primaryVariablesGlobal(transportEqIdx).size();
- TransportSolutionType transportedQuantity(problem.variables().primaryVariablesGlobal(transportEqIdx));
- TransportSolutionType transValueOldIter(problem.variables().primaryVariablesGlobal(transportEqIdx));
- TransportSolutionType transValueHelp(transSize);
- TransportSolutionType transValueDiff(transSize);
- TransportSolutionType updateOldIter(transSize);
- TransportSolutionType updateHelp(transSize);
- TransportSolutionType updateDiff(transSize);
-
bool converg = false;
int iter = 0;
int iterTot = 0;
- updateOldIter = 0;
while (!converg)
{
iter++;
iterTot++;
- problem.pressureModel().pressure(false);
+ problem.pressureModel().update(false);
//calculate velocities
+ // TODO: remove calculateVelocity in impet
problem.pressureModel().calculateVelocity();
//calculate defect of transported quantity
problem.transportModel().update(t, dt, updateVec, true);
- if (iterFlag_)
- { // only needed if iteration has to be done
+ if (iterFlag_)// only needed if iteration has to be done
+ {
+ // the method is valid for any transported quantity.
+ int transSize = problem.transportModel().transportedQuantity().size();
+ TransportSolutionType transportedQuantity(problem.transportModel().transportedQuantity());
+ TransportSolutionType transValueOldIter(problem.transportModel().transportedQuantity());
+ TransportSolutionType transValueHelp(transSize);
+ TransportSolutionType transValueDiff(transSize);
+ TransportSolutionType updateOldIter(transSize);
+ TransportSolutionType updateHelp(transSize);
+ TransportSolutionType updateDiff(transSize);
+
+ updateOldIter = 0;
updateHelp = updateVec;
- transportedQuantity = problem.variables().primaryVariablesGlobal(transportEqIdx);
+ transportedQuantity = problem.transportModel().transportedQuantity();
transportedQuantity += (updateHelp *= (dt * cFLFactor_));
transportedQuantity *= omega_;
transValueHelp = transValueOldIter;
@@ -144,30 +145,31 @@ public:
updateDiff -= updateOldIter;
transValueOldIter = transportedQuantity;
updateOldIter = updateVec;
+
+ // break criteria for iteration loop
+ if (iterFlag_ == 2 && dt * updateDiff.two_norm() / transportedQuantity.two_norm() <= maxDefect_)
+ {
+ converg = true;
+ }
+ else if (iterFlag_ == 1 && iter > nIter_)
+ {
+ converg = true;
+ }
+
+ if (iterFlag_ == 2 && transportedQuantity.infinity_norm() > (1 + maxDefect_))
+ {
+ converg = false;
+ }
+ if (!converg && iter > nIter_)
+ {
+ converg = true;
+ std::cout << "Nonlinear loop in IMPET.update exceeded nIter = " << nIter_ << " iterations." << std::endl;
+ std::cout << transportedQuantity.infinity_norm() << std::endl;
+ }
}
- // break criteria for iteration loop
- if (iterFlag_ == 2 && dt * updateDiff.two_norm() / transportedQuantity.two_norm() <= maxDefect_)
- {
- converg = true;
- }
- else if (iterFlag_ == 1 && iter > nIter_)
- {
- converg = true;
- }
- else if (iterFlag_ == 0)
- {
- converg = true;
- }
- if (iterFlag_ == 2 && transportedQuantity.infinity_norm() > (1 + maxDefect_))
- {
- converg = false;
- }
- if (!converg && iter > nIter_)
- {
+ else // no iteration => always "converged"
converg = true;
- std::cout << "Nonlinear loop in IMPET.update exceeded nIter = " << nIter_ << " iterations." << std::endl;
- std::cout << transportedQuantity.infinity_norm() << std::endl;
- }
+
}
// outputs
if (iterFlag_ == 2)
diff --git a/dumux/decoupled/common/impetproblem.hh b/dumux/decoupled/common/impetproblem.hh
index ced54283944f4b2d946cabe00c84e883c0026d09..d154e00026e111f11f89b8818c3ca1a3325386d9 100644
--- a/dumux/decoupled/common/impetproblem.hh
+++ b/dumux/decoupled/common/impetproblem.hh
@@ -77,8 +77,8 @@ private:
};
enum
{
- adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid),
- transportEqIdx = Indices::transportEqIdx
+ wetting = 0, nonwetting = 1,
+ adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid)
};
typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
@@ -411,14 +411,14 @@ public:
void timeIntegration()
{
// allocate temporary vectors for the updates
- typedef TransportSolutionType Solution;
- Solution k1 = asImp_().variables().primaryVariablesGlobal(transportEqIdx);
+ TransportSolutionType updateVector
+ = transportModel().transportedQuantity();
Scalar t = timeManager().time();
Scalar dt = 1e100;
// obtain the first update and the time step size
- model().update(t, dt, k1);
+ model().update(t, dt, updateVector);
//make sure t_old + dt is not larger than tend
dt = std::min(dt, timeManager().episodeMaxTimeStepSize());
@@ -449,13 +449,7 @@ public:
timeManager().setTimeStepSize(dt);
// explicit Euler: Sat <- Sat + dt*N(Sat)
- int size = gridView_.size(0);
- Solution& newSol = asImp_().variables().primaryVariablesGlobal(transportEqIdx);
- for (int i = 0; i < size; i++)
- {
- newSol[i] += (k1[i] *= dt);
- transportModel().updateSaturationSolution(i, newSol[i][0]);
- }
+ transportModel().updateTransportedQuantity(updateVector);
}
/*!
@@ -744,7 +738,10 @@ public:
timeManager().serialize(res);
resultWriter().serialize(res);
- model().serialize(res);
+
+ // do the actual serialization process: write primary variables
+ res.template serializeEntities<0> (*pressModel_, gridView_);
+ res.template serializeEntities<0> (*transportModel_, gridView_);
res.serializeEnd();
}
@@ -766,7 +763,11 @@ public:
timeManager().deserialize(res);
resultWriter().deserialize(res);
- model().deserialize(res);
+
+ // do the actual serialization process: get primary variables
+ res.template deserializeEntities<0> (*pressModel_, gridView_);
+ res.template deserializeEntities<0> (*transportModel_, gridView_);
+ pressureModel().updateMaterialLaws();
res.deserializeEnd();
};
diff --git a/dumux/decoupled/common/variableclass.hh b/dumux/decoupled/common/variableclass.hh
index d1262c8b0a57d8e30df2fe12dd17475dec5c6e2c..76a57a3023eeee06a0dc9722655e2a794025e64c 100644
--- a/dumux/decoupled/common/variableclass.hh
+++ b/dumux/decoupled/common/variableclass.hh
@@ -72,6 +72,7 @@ private:
public:
typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;//!<type for vector of scalars
+ typedef typename GET_PROP_TYPE(TypeTag, TransportSolutionType) TransportSolutionType;
typedef typename std::vector <CellData> CellDataVector;
private:
@@ -80,8 +81,6 @@ private:
VertexMapper vertexMapper_;
const int codim_;
CellDataVector cellDataVector_;
-protected:
- ScalarSolutionType primaryVariablesVector_[numEq];
public:
//! Constructs a VariableClass object
@@ -93,50 +92,9 @@ public:
VariableClass(const GridView& gridView, int codim = 0) :
gridView_(gridView), elementMapper_(gridView), vertexMapper_(gridView), codim_(codim)
{
- for (int i = 0; i < numEq; i++)
- {
- primaryVariablesVector_[i].resize(gridView.size(codim));
- primaryVariablesVector_[i] = 0;
- }
cellDataVector_.resize(gridView.size(codim));
}
- // serialization methods
- //! Function needed for restart option.
- template<class Restarter>
- void serialize(Restarter &res)
- {
- res.template serializeEntities<0> (*this, gridView_);
- }
-
- //! Function needed for restart option.
- template<class Restarter>
- void deserialize(Restarter &res)
- {
- res.template deserializeEntities<0> (*this, gridView_);
- }
-
- //! Function needed for restart option.
- void serializeEntity(std::ostream &outstream, const Element &element)
- {
- int globalIdx = elementMapper_.map(element);
- outstream << primaryVariablesVector_[0][globalIdx][0];
- for (int i = 1; i < numEq; i++)
- {
- outstream <<" "<< primaryVariablesVector_[i][globalIdx][0];
- }
- }
-
- //! Function needed for restart option.
- void deserializeEntity(std::istream &instream, const Element &element)
- {
- int globalIdx = elementMapper_.map(element);
- instream >> primaryVariablesVector_[0][globalIdx][0];
- for (int i = 1; i < numEq; i++)
- {
- instream >> primaryVariablesVector_[i][globalIdx][0];
- }
- }
//! Resizes decoupled variable vectors
/*! Method that change the size of the vectors for h-adaptive simulations.
@@ -144,37 +102,12 @@ public:
*\param size Size of the current (refined and coarsened) grid
*/
void adaptVariableSize(int size)
- {
- for (int i = 0; i < numEq; i++)
- {
- primaryVariablesVector_[i].resize(size);
- primaryVariablesVector_[i] = 0;
- }
- cellDataVector_.resize(size);
- }
-
-// void communicatePressure()
-// {
-//#if HAVE_MPI
-// ElementHandleAssign<typename ScalarSolutionType::block_type, ScalarSolutionType, ElementMapper> elementHandle(primaryVariables_, elementMapper_);
-// gridView_.communicate(elementHandle,
-// Dune::InteriorBorder_All_Interface,
-// Dune::ForwardCommunication);
-//#endif
-// }
-
- //! Return pressure vector
- const ScalarSolutionType& primaryVariablesGlobal(int eqIdx) const
{
- return primaryVariablesVector_[eqIdx];
- }
-
- ScalarSolutionType& primaryVariablesGlobal(int eqIdx)
- {
- return primaryVariablesVector_[eqIdx];
- }
+ DUNE_THROW(Dune::NotImplemented,"TODO: Primary Variables Vector have to be resized!!!");
+ cellDataVector_.resize(size);
+ }
- //! Return saturation vector
+ //! Return the vector holding all cell data
CellDataVector& cellDataGlobal()
{
return cellDataVector_;
@@ -185,7 +118,7 @@ public:
return cellDataVector_;
}
- //! Return saturation vector
+ //! Return the cell data of a specific cell
CellData& cellData(int idx)
{
return cellDataVector_[idx];