diff --git a/dumux/implicit/adaptive/gridadapt.hh b/dumux/implicit/adaptive/gridadapt.hh
index 2c9f327089401698a42967e7375d0ab50cfb28d5..cee91b80d2702e55b824825ccd4a664cb7b592af 100644
--- a/dumux/implicit/adaptive/gridadapt.hh
+++ b/dumux/implicit/adaptive/gridadapt.hh
@@ -77,14 +77,6 @@ public:
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
-// {
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 +85,6 @@ public:
{
DUNE_THROW(Dune::InvalidStateException, "Coarsening the level 0 entities is not possible! Choose MinLevel >= 0");
}
-// }
}
/*!
@@ -149,9 +140,77 @@ public:
*/
void adaptGrid()
{
- adaptGrid(adaptionIndicator_) ;
+ if(levelMax_ > levelMin_)
+ adaptGrid(adaptionIndicator_) ;
}
+ /*!
+ * @brief Returns true if grid cells have been marked for adaption
+ */
+ bool wasAdapted()
+ {
+ int sumMarked = problem_.grid().comm().sum(marked_);
+ int sumCoarsened = problem_.grid().comm().sum(coarsened_);
+
+ return (sumMarked != 0 || sumCoarsened != 0);
+ }
+
+ const bool wasAdapted() const
+ {
+ int sumMarked = problem_.grid().comm().sum(marked_);
+ int sumCoarsened = problem_.grid().comm().sum(coarsened_);
+
+ return (sumMarked != 0 || sumCoarsened != 0);
+ }
+
+ /*!
+ * Sets minimum and maximum refinement levels
+ *
+ * @param levMin minimum level for coarsening
+ * @param levMax maximum level for refinement
+ */
+ void setLevels(int levMin, int levMax)
+ {
+ if (levMin < 0)
+ DUNE_THROW(Dune::InvalidStateException, "Coarsening the level 0 entities is not possible!");
+ levelMin_ = levMin;
+ levelMax_ = levMax;
+ }
+
+ /*!
+ * @brief Returns maximum refinement level
+ *
+ * The value is the assign maximum possible level,
+ * not the actual maximum level of the grid.
+ * @return levelMax_ maximum level for refinement
+ */
+ const int getMaxLevel() const
+ {
+ return levelMax_;
+ }
+ /*!
+ * @brief Returns minimum refinement level
+ *
+ * The value is the assign minimum possible level,
+ * not the actual minimum level of the grid.
+ * @return levelMin_ minimum level for coarsening
+ */
+ const int getMinLevel() const
+ {
+ return levelMin_;
+ }
+
+ AdaptionIndicator& adaptionIndicator()
+ {
+ return adaptionIndicator_;
+ }
+
+ AdaptionIndicator& adaptionIndicator() const
+ {
+ return adaptionIndicator_;
+ }
+
+private:
/*!
* @brief Method to adapt the grid with individual indicator vector
*
@@ -257,73 +316,6 @@ public:
}
}
- /*!
- * @brief Returns true if grid cells have been marked for adaption
- */
- bool wasAdapted()
- {
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
-
- return (sumMarked != 0 || sumCoarsened != 0);
- }
-
- const bool wasAdapted() const
- {
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
-
- return (sumMarked != 0 || sumCoarsened != 0);
- }
-
- /*!
- * Sets minimum and maximum refinement levels
- *
- * @param levMin minimum level for coarsening
- * @param levMax maximum level for refinement
- */
- void setLevels(int levMin, int levMax)
- {
- if (levMin < 0)
- DUNE_THROW(Dune::InvalidStateException, "Coarsening the level 0 entities is not possible!");
- levelMin_ = levMin;
- levelMax_ = levMax;
- }
-
- /*!
- * @brief Returns maximum refinement level
- *
- * The value is the assign maximum possible level,
- * not the actual maximum level of the grid.
- * @return levelMax_ maximum level for refinement
- */
- const int getMaxLevel() const
- {
- return levelMax_;
- }
- /*!
- * @brief Returns minimum refinement level
- *
- * The value is the assign minimum possible level,
- * not the actual minimum level of the grid.
- * @return levelMin_ minimum level for coarsening
- */
- const int getMinLevel() const
- {
- return levelMin_;
- }
-
- AdaptionIndicator& adaptionIndicator()
- {
- return adaptionIndicator_;
- }
-
- AdaptionIndicator& adaptionIndicator() const
- {
- return adaptionIndicator_;
- }
-
-private:
/*!
* @brief Method ensuring the refinement ratio of 2:1
*
diff --git a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
index 1cf2118c1394f56600cb95e71cbf900cb9e8f5fa..300c7708023e36d3d50bb3091c904f073a7a656f 100644
--- a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
+++ b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
@@ -66,7 +66,6 @@ private:
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;
@@ -123,7 +122,7 @@ public:
//get your map entry
AdaptedValues &adaptedValues = adaptionMap_[element];
- // put your value in the map
+ // put values in the map
if (element.isLeaf())
{
FVElementGeometry fvGeometry;
@@ -157,9 +156,17 @@ public:
//Average in father
if (element.level() > 0)
{
+ if(!element.hasFather())
+ DUNE_THROW(Dune::InvalidStateException, "Element on level > 0 has no father element!");
+
AdaptedValues& adaptedValuesFather = adaptionMap_[element.father()];
- adaptedValuesFather.count += 1;
- storeAdaptionValues(adaptedValues, adaptedValuesFather);
+ //For some grids the father element is identical to the son element.
+ //For that case averaging is not necessary.
+ if(&adaptedValues != &adaptedValuesFather)
+ {
+ adaptedValuesFather.count += 1;
+ storeAdaptionValues(adaptedValues, adaptedValuesFather);
+ }
}
if(isBox && !element.isLeaf())
@@ -180,40 +187,6 @@ public:
}
}
- 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)
*
@@ -249,7 +222,7 @@ public:
if (element.partitionType() == Dune::GhostEntity)
continue;
- if (!element.isNew())
+ if (!element.isNew() || element.level() == 0)
{
//entry is in map, write in leaf
if (element.isLeaf())
@@ -263,7 +236,7 @@ public:
{
// get index
auto subEntity = element.template subEntity <dofCodim>(scvIdx);
- int dofIdx = problem.model().dofMapper().subIndex(element,scvIdx,dofCodim);
+ int dofIdx = problem.model().dofMapper().index(subEntity);
this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdx],scvIdx);
@@ -299,19 +272,16 @@ public:
else
{
// value is not in map, interpolate from father element
- if (element.level() > 0 && element.hasFather())
+ if (element.hasFather())
{
auto eFather = element.father();
+ while(eFather.isNew() && eFather.level() > 0) eFather = eFather.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)
@@ -427,13 +397,16 @@ public:
}
}
}
+ else
+ {
+ DUNE_THROW(Dune::InvalidStateException, "Element is new but has no father element!");
+ }
}
}
}
- std::cout << "End of reconstructon" << std::endl;
if(isBox)
{
for(int dofIdx = 0; dofIdx < problem.model().numDofs(); dofIdx++)
@@ -442,7 +415,6 @@ public:
}
}
-
// reset entries in restrictionmap
adaptionMap_.resize( typename PersistentContainer::Value() );
adaptionMap_.shrinkToFit();