diff --git a/dumux/implicit/adaptive/adaptionhelper.hh b/dumux/implicit/adaptive/adaptionhelper.hh
index e88021c201076db895420931d1d93f0edbc1ae75..36ddbb7eec79deb3f05c5bb718fdd4b24df7f679 100644
--- a/dumux/implicit/adaptive/adaptionhelper.hh
+++ b/dumux/implicit/adaptive/adaptionhelper.hh
@@ -37,9 +37,7 @@ namespace Properties
{
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(ImplicitIsBox);
-NEW_PROP_TAG(PrimaryVariables);
NEW_PROP_TAG(Problem);
-NEW_PROP_TAG(Scalar);
}
/*!
@@ -51,47 +49,20 @@ class ImplicitAdaptionHelper
private:
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::Grid Grid;
enum {
// Grid and world dimension
dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
+ dimWorld = GridView::dimensionworld,
};
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
- typedef typename GridView::Grid Grid;
- typedef typename Grid::LevelGridView LevelGridView;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Traits::template Codim<dofCodim>::Entity DofEntity;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename GridView::ctype CoordScalar;
- typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
-
- typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
- typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
-
-protected:
- const GridView gridView_;
- const Grid& grid_;
-
- struct AdaptedValues
- {
- PrimaryVariables u;
- int count;
- AdaptedValues()
- {
- count = 0;
- }
- };
-
- typedef Dune::PersistentContainer<Grid, AdaptedValues> PersistentContainer;
- PersistentContainer adaptionMap_;
-
public:
//! Constructs an adaptive helper object
/**
@@ -100,11 +71,9 @@ public:
*
* @param gridView a DUNE gridview object corresponding to diffusion and transport equation
*/
- ImplicitAdaptionHelper(const GridView& gridView) :
- gridView_(gridView), grid_(gridView.grid()), adaptionMap_(grid_, dofCodim)
+ ImplicitAdaptionHelper(Problem& problem)
{}
-
/*!
* Store primary variables
*
@@ -118,60 +87,12 @@ public:
*/
void storePrimVars(Problem& problem)
{
- adaptionMap_.resize();
-
- // loop over all levels of the grid
- for (int level = grid_.maxLevel(); level >= 0; level--)
- {
- //get grid view on level grid
- LevelGridView levelView = grid_.levelGridView(level);
-
- if(!isBox)
- {
- for (const auto& element : elements(levelView))
- {
- //get your map entry
- AdaptedValues &adaptedValues = adaptionMap_[element];
-
- // put your value in the map
- if (element.isLeaf())
- {
- // get index
- int indexI = this->elementIndex(problem, element);
-
- storeAdaptionValues(adaptedValues, problem.model().curSol()[indexI]);
-
- adaptedValues.count = 1;
- }
- //Average in father
- if (element.level() > 0)
- {
- AdaptedValues& adaptedValuesFather = adaptionMap_[element.father()];
- adaptedValuesFather.count += 1;
- storeAdaptionValues(adaptedValues, adaptedValuesFather);
- }
-
- }
- }
- else
- {
- for (const auto& entity : entities(levelView, Dune::Codim<dofCodim>()))
- {
- //get your map entry
- AdaptedValues &adaptedValues = adaptionMap_[entity];
-
- // put your value in the map
- int indexI = this->dofIndex(problem, entity);
-
- storeAdaptionValues(adaptedValues, problem.model().curSol()[indexI]);
-
- adaptedValues.count = 1;
-
- }
- }
- }
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide "
+ "a storePrimVars method.");
}
+
/*!
* Reconstruct missing primary variables (where elements are created/deleted)
*
@@ -187,207 +108,21 @@ public:
*/
void reconstructPrimVars(Problem& problem)
{
- adaptionMap_.resize();
-
- for (int level = 0; level <= grid_.maxLevel(); level++)
- {
- LevelGridView levelView = grid_.levelGridView(level);
-
- for (const auto& element : elements(levelView))
- {
- // only treat non-ghosts, ghost data is communicated afterwards
- if (element.partitionType() == Dune::GhostEntity)
- continue;
-
- if (!element.isNew())
- {
- //entry is in map, write in leaf
- if (element.isLeaf())
- {
- if(!isBox)
- {
- AdaptedValues &adaptedValues = adaptionMap_[element];
- int newIdxI = this->elementIndex(problem, element);
-
- setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
- }
- else
- {
- unsigned int numSubEntities = element.subEntities(dofCodim);
-
- for(unsigned int i = 0; i < numSubEntities; i++)
- {
- auto subEntity = element.template subEntity <dofCodim>(i);
- AdaptedValues &adaptedValues = adaptionMap_[subEntity];
- int newIdxI = this->dofIndex(problem, subEntity);
-
- setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
-
- }
- }
- }
- }
- else
- {
- // value is not in map, interpolate from father element
- if (element.level() > 0 && element.hasFather())
- {
- auto epFather = element.father();
-
- if(!isBox)
- {
- // create new entry: reconstruct from adaptionMap_[*father] to a new
- // adaptionMap_[*son]
- reconstructAdaptionValues(adaptionMap_, epFather, element, problem);
-
- // access new son
- AdaptedValues& adaptedValues = adaptionMap_[element];
- adaptedValues.count = 1;
-
- // if we are on leaf, store reconstructed values of son in CellData object
- if (element.isLeaf())
- {
- // acess new CellData object
- int newIdxI = this->elementIndex(problem, element);
-
- setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
- }
- }
- else
- {
- unsigned int numSubEntities= element.subEntities(dofCodim);
- const auto geometryI = element.geometry();
-
- for(unsigned int i = 0; i < numSubEntities; i++)
- {
- auto subEntity = element.template subEntity <dofCodim>(i);
- AdaptedValues &adaptedValues = adaptionMap_[subEntity];
-
- if(adaptedValues.count == 0){
- LocalPosition dofCenterPos = geometryI.local(subEntity.geometry().center());
- const LocalFiniteElementCache feCache;
- Dune::GeometryType geomType = epFather.geometry().type();
-
- const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
- std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
- localFiniteElement.localBasis().evaluateFunction(dofCenterPos, shapeVal);
- PrimaryVariables u(0);
- for (int j = 0; j < shapeVal.size(); ++j)
- {
- AdaptedValues & adaptedValuesFather = adaptionMap_[epFather.template subEntity <dofCodim>(j)];
- u.axpy(shapeVal[j], adaptedValuesFather.u);
- }
-
- adaptedValues.u = u;
- adaptedValues.count = 1;
- }
-
- if (element.isLeaf())
- {
- int newIdxI = this->dofIndex(problem, subEntity);
- setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]);
- }
-
- }
- }
- }
-
- }
- }
-
- }
- // reset entries in restrictionmap
- adaptionMap_.resize( typename PersistentContainer::Value() );
- adaptionMap_.shrinkToFit();
- adaptionMap_.fill( typename PersistentContainer::Value() );
-
-//#if HAVE_MPI
-// // communicate ghost data
-// typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
-// typedef typename SolutionTypes::ElementMapper ElementMapper;
-// typedef VectorExchange<ElementMapper, std::vector<CellData> > DataHandle;
-// DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
-// problem.gridView().template communicate<DataHandle>(dataHandle,
-// Dune::InteriorBorder_All_Interface,
-// Dune::ForwardCommunication);
-//#endif
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide "
+ "a reconstructPrimVars method.");
}
- //! Stores values to be adapted in an adaptedValues container
- /**
- * Stores values to be adapted from the current CellData objects into
- * the adaption container in order to be mapped on a new grid.
- *
- * \param adaptedValues Container for model-specific values to be adapted
- * \param u The variables to be stored
- */
- static void storeAdaptionValues(AdaptedValues& adaptedValues, const PrimaryVariables& u)
- {
- adaptedValues.u = u;
- }
- //! Stores sons entries into father element for averaging
- /**
- * Sum up the adaptedValues (sons values) into father element. We store from leaf
- * upwards, so sons are stored first, then cells on the next leaf (=fathers)
- * can be averaged.
- *
- * \param adaptedValues Container for model-specific values to be adapted
- * \param adaptedValuesFather Values to be adapted of father cell
- */
- static void storeAdaptionValues(AdaptedValues& adaptedValues,
- AdaptedValues& adaptedValuesFather)
- {
- if(!isBox)
- {
- adaptedValuesFather.u += adaptedValues.u;
- adaptedValuesFather.u /= adaptedValues.count;
- }
- else
- {
- adaptedValuesFather.u = adaptedValues.u;
- }
- }
- //! Set adapted values in CellData
- /**
- * This methods stores reconstructed values into the cellData object, by
- * this setting a newly mapped solution to the storage container of the
- * sequential models.
- *
- * \param adaptedValues Container for model-specific values to be adapted
- * \param u The variables to be stored
- */
- static void setAdaptionValues(AdaptedValues& adaptedValues, PrimaryVariables& u)
- {
- PrimaryVariables uNew = adaptedValues.u;
- uNew /= adaptedValues.count;
-
- u = uNew;
- }
+protected:
- //! Reconstructs sons entries from data of father cell
- /**
- * Reconstructs a new solution from a father cell into a newly
- * generated son cell. New cell is stored into the global
- * adaptionMap.
- *
- * \param adaptionMap Global map storing all values to be adapted
- * \param father Entity Pointer to the father cell
- * \param son Entity Pointer to the newly created son cell
- * \param problem The problem
- */
- static void reconstructAdaptionValues(Dune::PersistentContainer<Grid, AdaptedValues>& adaptionMap,
- const Element& father, const Element& son, const Problem& problem)
+ int dofIndex(const Problem& problem, const DofEntity& entity) const
{
- AdaptedValues& adaptedValues = adaptionMap[son];
- AdaptedValues& adaptedValuesFather = adaptionMap[father];
-
- adaptedValues.u = adaptedValuesFather.u;
- adaptedValues.u /= adaptedValuesFather.count;
+ return problem.model().dofMapper().index(entity);
}
- int dofIndex(const Problem& problem, const DofEntity& entity) const
+ int dofIndex(const Problem& problem, const Element& element, const int scvIdx) const
{
- return problem.model().dofMapper().index(entity);
+ return problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
}
int elementIndex(const Problem& problem, const Element& element) const
diff --git a/dumux/implicit/adaptive/gridadapt.hh b/dumux/implicit/adaptive/gridadapt.hh
index ccc221d6609df4b478f31c471fdd4bf421770459..2c9f327089401698a42967e7375d0ab50cfb28d5 100644
--- a/dumux/implicit/adaptive/gridadapt.hh
+++ b/dumux/implicit/adaptive/gridadapt.hh
@@ -72,19 +72,19 @@ public:
*/
ImplicitGridAdapt (Problem& problem)
: problem_(problem),
- adaptionHelper_(problem.gridView()),
+ adaptionHelper_(problem),
adaptionIndicator_(problem),
marked_(0),
coarsened_(0)
{
- if(isBox)
- {
- DUNE_THROW(Dune::NotImplemented,
- "Grid adaption is not yet mass conservative for Box method! "
- << "Use cell-centered scheme instead!");
- }
- else
- {
+// if(isBox)
+// {
+// DUNE_THROW(Dune::NotImplemented,
+// "Grid adaption is not yet mass conservative for Box method! "
+// << "Use cell-centered scheme instead!");
+// }
+// else
+// {
levelMin_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MinLevel);
levelMax_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MaxLevel);
adaptionInterval_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, AdaptionInterval);
@@ -93,7 +93,7 @@ public:
{
DUNE_THROW(Dune::InvalidStateException, "Coarsening the level 0 entities is not possible! Choose MinLevel >= 0");
}
- }
+// }
}
/*!
diff --git a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1cf2118c1394f56600cb95e71cbf900cb9e8f5fa
--- /dev/null
+++ b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
@@ -0,0 +1,539 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * 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_TWOP_ADAPTIONHELPER_HH
+#define DUMUX_TWOP_ADAPTIONHELPER_HH
+
+#include <dumux/implicit/adaptive/adaptionhelper.hh>
+
+namespace Dumux {
+
+namespace Properties
+{
+NEW_PROP_TAG(PrimaryVariables);
+NEW_PROP_TAG(Scalar);
+}
+
+/*!
+ * \brief Base class holding the variables for implicit models.
+ */
+template<class TypeTag>
+class TwoPAdaptionHelper : public ImplicitAdaptionHelper<TypeTag>
+{
+private:
+ typedef ImplicitAdaptionHelper<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ //indices
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pwsn = Indices::pwsn,
+ pnsw = Indices::pnsw,
+ formulation = GET_PROP_VALUE(TypeTag, Formulation),
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+ typedef typename GridView::Grid Grid;
+ typedef typename Grid::LevelGridView LevelGridView;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Traits::template Codim<dofCodim>::Entity DofEntity;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef typename GridView::ctype CoordScalar;
+ typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+
+ typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
+ typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
+
+ struct AdaptedValues
+ {
+ std::vector<PrimaryVariables> u;
+ int count;
+ PrimaryVariables associatedMass;
+ AdaptedValues(): associatedMass(0)
+ {
+ count = 0;
+ }
+ };
+
+ typedef Dune::PersistentContainer<Grid, AdaptedValues> PersistentContainer;
+ PersistentContainer adaptionMap_;
+
+public:
+ //! Constructs an adaption helper object
+ /**
+ * @param gridView a DUNE gridview object
+ */
+ TwoPAdaptionHelper(Problem& problem) : ParentType(problem), adaptionMap_(problem.grid(), 0)
+ {}
+
+ /*!
+ * Store primary variables
+ *
+ * To reconstruct the solution in father elements, problem properties might
+ * need to be accessed.
+ * From upper level on downwards, the old solution is stored into an container
+ * object, before the grid is adapted. Father elements hold averaged information
+ * from the son cells for the case of the sons being coarsened.
+ *
+ * @param problem The current problem
+ */
+ void storePrimVars(Problem& problem)
+ {
+ adaptionMap_.resize();
+
+ // loop over all levels of the grid
+ for (int level = problem.grid().maxLevel(); level >= 0; level--)
+ {
+ //get grid view on level grid
+ LevelGridView levelView = problem.grid().levelGridView(level);
+
+ for (const auto& element : elements(levelView))
+ {
+ //get your map entry
+ AdaptedValues &adaptedValues = adaptionMap_[element];
+
+ // put your value in the map
+ if (element.isLeaf())
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem.gridView(), element);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ // get index
+ auto subEntity = element.template subEntity <dofCodim>(scvIdx);
+ int dofIdx = this->dofIndex(problem, subEntity);
+
+ adaptedValues.u.push_back(problem.model().curSol()[dofIdx]);
+
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u[scvIdx],
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ false);
+
+ Scalar volume = fvGeometry.subContVol[scvIdx].volume;
+
+ adaptedValues.associatedMass[nPhaseIdx] += volume*volVars.density(nPhaseIdx)
+ * volVars.porosity() * volVars.saturation(nPhaseIdx);
+ adaptedValues.associatedMass[wPhaseIdx] += volume*volVars.density(wPhaseIdx)
+ * volVars.porosity() * volVars.saturation(wPhaseIdx);
+ }
+ adaptedValues.count = 1;
+ }
+ //Average in father
+ if (element.level() > 0)
+ {
+ AdaptedValues& adaptedValuesFather = adaptionMap_[element.father()];
+ adaptedValuesFather.count += 1;
+ storeAdaptionValues(adaptedValues, adaptedValuesFather);
+ }
+
+ if(isBox && !element.isLeaf())
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem.gridView(), element);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ auto subEntity = element.template subEntity <dofCodim>(scvIdx);
+ int dofIdx = this->dofIndex(problem, subEntity);
+
+ adaptedValues.u.push_back(problem.model().curSol()[dofIdx]);
+ }
+ }
+
+ }
+ }
+ }
+
+ void dataOutput(Problem& problem)
+ {
+ // loop over all levels of the grid
+ for (int level = problem.grid().maxLevel(); level >= 0; level--)
+ {
+ //get grid view on level grid
+ LevelGridView levelView = problem.grid().levelGridView(level);
+
+ for (const auto& element : elements(levelView))
+ {
+ AdaptedValues &adaptedValues = adaptionMap_[element];
+ if (element.level() > 0 && !element.isNew())
+ {
+ AdaptedValues& adaptedValuesFather = adaptionMap_[element.father()];
+ for(int i=0; i<adaptedValues.u.size(); i++)
+ {
+ auto subEntity = element.template subEntity <dofCodim>(i);
+ int dofIdx = this->dofIndex(problem, subEntity);
+ std::cout << "SonValues:" << std::endl;
+ std::cout << "u: " << dofIdx <<"___" << adaptedValues.u[i] << std::endl;
+ }
+ for(int i=0; i<adaptedValuesFather.u.size(); i++)
+ {
+ auto subEntity = element.father().template subEntity <dofCodim>(i);
+ int dofIdx = this->dofIndex(problem, subEntity);
+ std::cout << "FatherValues:" << std::endl;
+ std::cout << "u: " << dofIdx <<"___" << adaptedValuesFather.u[i] << std::endl;
+ }
+ }
+
+ }
+ }
+ }
+
+ /*!
+ * Reconstruct missing primary variables (where elements are created/deleted)
+ *
+ * To reconstruct the solution in father elements, problem properties might
+ * need to be accessed.
+ * Starting from the lowest level, the old solution is mapped on the new grid:
+ * Where coarsened, new cells get information from old father element.
+ * Where refined, a new solution is reconstructed from the old father cell,
+ * and then a new son is created. That is then stored into the general data
+ * structure (CellData).
+ *
+ * @param problem The current problem
+ */
+ void reconstructPrimVars(Problem& problem)
+ {
+ adaptionMap_.resize();
+ std::vector<Scalar> massCoeff;
+ std::vector<Scalar> associatedMass;
+
+ if(isBox)
+ {
+ massCoeff.resize(problem.model().numDofs(),0.0);
+ associatedMass.resize(problem.model().numDofs(),0.0);
+ }
+
+ for (int level = 0; level <= problem.grid().maxLevel(); level++)
+ {
+ LevelGridView levelView = problem.grid().levelGridView(level);
+
+ for (const auto& element : elements(levelView))
+ {
+ // only treat non-ghosts, ghost data is communicated afterwards
+ if (element.partitionType() == Dune::GhostEntity)
+ continue;
+
+ if (!element.isNew())
+ {
+ //entry is in map, write in leaf
+ if (element.isLeaf())
+ {
+ AdaptedValues &adaptedValues = adaptionMap_[element];
+
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem.gridView(), element);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ // get index
+ auto subEntity = element.template subEntity <dofCodim>(scvIdx);
+ int dofIdx = problem.model().dofMapper().subIndex(element,scvIdx,dofCodim);
+
+ this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdx],scvIdx);
+
+ if(isBox)
+ {
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u[scvIdx],
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ false);
+
+ Scalar volume = fvGeometry.subContVol[scvIdx].volume;
+ Scalar volumeElement = fvGeometry.elementVolume;
+
+ if (int(formulation) == pwsn)
+ {
+ massCoeff[dofIdx] += volume * volVars.density(nPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdx] += volume/volumeElement*adaptedValues.associatedMass[nPhaseIdx];
+ }
+ else if (int(formulation) == pnsw)
+ {
+ massCoeff[dofIdx] += volume * volVars.density(wPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdx] += volume/volumeElement*adaptedValues.associatedMass[wPhaseIdx];
+ }
+ }
+
+ }
+
+ }
+ }
+ else
+ {
+ // value is not in map, interpolate from father element
+ if (element.level() > 0 && element.hasFather())
+ {
+ auto eFather = element.father();
+ FVElementGeometry fvGeometryFather;
+ fvGeometryFather.update(problem.gridView(), eFather);
+ Scalar massFather = 0.0;
+
+ if(!isBox)
+ {
+ // create new entry: reconstruct from adaptionMap_[*father] to a new
+ // adaptionMap_[*son]
+ //this->reconstructAdaptionValues(this->adaptionMap_, eFather, element, problem);
+
+ AdaptedValues& adaptedValuesFather = adaptionMap_[eFather];
+
+ if (int(formulation) == pwsn)
+ {
+ massFather = adaptedValuesFather.associatedMass[nPhaseIdx];
+ }
+ else if (int(formulation) == pnsw)
+ {
+ massFather = adaptedValuesFather.associatedMass[wPhaseIdx];
+ }
+
+ // access new son
+ AdaptedValues& adaptedValues = adaptionMap_[element];
+ adaptedValues.count = 1;
+
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem.gridView(), element);
+
+ adaptedValues.u.push_back(adaptedValuesFather.u[0]);
+
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u[0],
+ problem,
+ element,
+ fvGeometry,
+ /*scvIdx=*/0,
+ false);
+
+ Scalar volume = fvGeometry.subContVol[0].volume;
+ Scalar massCoeffSon = 0.0;
+ if (int(formulation) == pwsn)
+ {
+ massCoeffSon = volume * volVars.density(nPhaseIdx) * volVars.porosity();
+ }
+ else if (int(formulation) == pnsw)
+ {
+ massCoeffSon = volume * volVars.density(wPhaseIdx) * volVars.porosity();
+ }
+ Scalar volumeFather = eFather.geometry().volume();
+ adaptedValues.u[0][saturationIdx] = volume/volumeFather*massFather/massCoeffSon;
+
+ // if we are on leaf, store reconstructed values of son in CellData object
+ if (element.isLeaf())
+ {
+ // acess new CellData object
+ int newIdxI = this->elementIndex(problem, element);
+
+ this->setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI],0);
+ }
+ }
+ else
+ {
+ AdaptedValues& adaptedValuesFather = adaptionMap_[eFather];
+ // access new son
+ AdaptedValues& adaptedValues = adaptionMap_[element];
+ adaptedValues.u.clear();
+ adaptedValues.count = 1;
+
+ const auto geometryI = element.geometry();
+
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem.gridView(), element);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ auto subEntity = element.template subEntity <dofCodim>(scvIdx);
+
+ LocalPosition dofCenterPos = geometryI.local(subEntity.geometry().center());
+ const LocalFiniteElementCache feCache;
+ Dune::GeometryType geomType = eFather.geometry().type();
+
+ const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
+ std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
+ localFiniteElement.localBasis().evaluateFunction(dofCenterPos, shapeVal);
+ PrimaryVariables u(0);
+ for (int j = 0; j < shapeVal.size(); ++j)
+ {
+ u.axpy(shapeVal[j], adaptedValuesFather.u[j]);
+ }
+
+ adaptedValues.u.push_back(u);
+ adaptedValues.count = 1;
+
+ if (element.isLeaf())
+ {
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u[scvIdx],
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ false);
+
+ Scalar volume = fvGeometry.subContVol[scvIdx].volume;
+ Scalar volumeFather = eFather.geometry().volume();
+
+ //int dofIdx = this->dofIndex(problem, subEntity);
+ int dofIdx = problem.model().dofMapper().subIndex(element,scvIdx,dofCodim);
+ if (int(formulation) == pwsn)
+ {
+ massCoeff[dofIdx] += volume * volVars.density(nPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdx] += volume/volumeFather*adaptedValuesFather.associatedMass[nPhaseIdx];
+ }
+ else if (int(formulation) == pnsw)
+ {
+ massCoeff[dofIdx] += volume * volVars.density(wPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdx] += volume/volumeFather*adaptedValuesFather.associatedMass[wPhaseIdx];
+ }
+
+ this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdx],scvIdx);
+ }
+
+ }
+ }
+ }
+
+ }
+ }
+
+ }
+
+ std::cout << "End of reconstructon" << std::endl;
+ if(isBox)
+ {
+ for(int dofIdx = 0; dofIdx < problem.model().numDofs(); dofIdx++)
+ {
+ problem.model().curSol()[dofIdx][saturationIdx] = associatedMass[dofIdx]/massCoeff[dofIdx];
+ }
+ }
+
+
+ // reset entries in restrictionmap
+ adaptionMap_.resize( typename PersistentContainer::Value() );
+ adaptionMap_.shrinkToFit();
+ adaptionMap_.fill( typename PersistentContainer::Value() );
+
+//#if HAVE_MPI
+// // communicate ghost data
+// typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+// typedef typename SolutionTypes::ElementMapper ElementMapper;
+// typedef VectorExchange<ElementMapper, std::vector<CellData> > DataHandle;
+// DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
+// problem.gridView().template communicate<DataHandle>(dataHandle,
+// Dune::InteriorBorder_All_Interface,
+// Dune::ForwardCommunication);
+//#endif
+ }
+
+ //! Stores sons entries into father element for averaging
+ /**
+ * Sum up the adaptedValues (sons values) into father element. We store from leaf
+ * upwards, so sons are stored first, then cells on the next leaf (=fathers)
+ * can be averaged.
+ *
+ * \param adaptedValues Container for model-specific values to be adapted
+ * \param adaptedValuesFather Values to be adapted of father cell
+ */
+ static void storeAdaptionValues(AdaptedValues& adaptedValues,
+ AdaptedValues& adaptedValuesFather)
+ {
+ if(!isBox)
+ {
+ if(adaptedValuesFather.u.size() == 0) adaptedValuesFather.u.resize(1);
+
+ adaptedValuesFather.u[0] += adaptedValues.u[0];
+ adaptedValuesFather.u[0] /= adaptedValues.count;
+ adaptedValuesFather.associatedMass += adaptedValues.associatedMass;
+ }
+ else
+ {
+ //adaptedValuesFather.u = adaptedValues.u;
+ adaptedValuesFather.associatedMass += adaptedValues.associatedMass;
+ }
+ }
+ //! Set adapted values in CellData
+ /**
+ * This methods stores reconstructed values into the cellData object, by
+ * this setting a newly mapped solution to the storage container of the
+ * sequential models.
+ *
+ * \param adaptedValues Container for model-specific values to be adapted
+ * \param u The variables to be stored
+ */
+ static void setAdaptionValues(AdaptedValues& adaptedValues, PrimaryVariables& u, int scvIdx)
+ {
+ PrimaryVariables uNew(0);
+ if(!isBox)
+ {
+ uNew = adaptedValues.u[scvIdx];
+ uNew /= adaptedValues.count;
+ }
+ else
+ {
+ uNew = adaptedValues.u[scvIdx];
+ }
+
+ u = uNew;
+ }
+
+ //! Reconstructs sons entries from data of father cell
+ /**
+ * Reconstructs a new solution from a father cell into a newly
+ * generated son cell. New cell is stored into the global
+ * adaptionMap.
+ *
+ * \param adaptionMap Global map storing all values to be adapted
+ * \param father Entity Pointer to the father cell
+ * \param son Entity Pointer to the newly created son cell
+ * \param problem The problem
+ */
+ static void reconstructAdaptionValues(Dune::PersistentContainer<Grid, AdaptedValues>& adaptionMap,
+ const Element& father, const Element& son, const Problem& problem)
+ {
+ AdaptedValues& adaptedValues = adaptionMap[son];
+ AdaptedValues& adaptedValuesFather = adaptionMap[father];
+
+ adaptedValues.u = adaptedValuesFather.u;
+ adaptedValues.u /= adaptedValuesFather.count;
+ }
+
+
+
+};
+}
+#endif