diff --git a/dumux/implicit/adaptive/adaptionhelper.hh b/dumux/implicit/adaptive/adaptionhelper.hh
index 36ddbb7eec79deb3f05c5bb718fdd4b24df7f679..21eba4e4772484066df388c6ccbaa650e5bbfa4c 100644
--- a/dumux/implicit/adaptive/adaptionhelper.hh
+++ b/dumux/implicit/adaptive/adaptionhelper.hh
@@ -47,21 +47,7 @@ template<class TypeTag>
class ImplicitAdaptionHelper
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::Grid Grid;
-
- enum {
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- };
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
-
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::Traits::template Codim<dofCodim>::Entity DofEntity;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
public:
//! Constructs an adaptive helper object
@@ -113,23 +99,8 @@ public:
"a reconstructPrimVars method.");
}
-protected:
-
- int dofIndex(const Problem& problem, const DofEntity& entity) const
- {
- return problem.model().dofMapper().index(entity);
- }
-
- int dofIndex(const Problem& problem, const Element& element, const int scvIdx) const
- {
- return problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
- }
+};
- int elementIndex(const Problem& problem, const Element& element) const
- {
- return problem.elementMapper().index(element);
- }
+} // end namespace Dumux
-};
-}
#endif
diff --git a/dumux/implicit/adaptive/gridadapt.hh b/dumux/implicit/adaptive/gridadapt.hh
index 015369d71731d6ed1a5ecda031f4fb5942c6d83b..17c1af7ff4e54df10a20ca28f416aae1374e0fa1 100644
--- a/dumux/implicit/adaptive/gridadapt.hh
+++ b/dumux/implicit/adaptive/gridadapt.hh
@@ -74,8 +74,8 @@ public:
: problem_(problem),
adaptionHelper_(problem),
adaptionIndicator_(problem),
- marked_(0),
- coarsened_(0)
+ sumRefine_(0),
+ sumCoarsen_(0)
{
levelMin_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MinLevel);
levelMax_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MaxLevel);
@@ -149,10 +149,7 @@ public:
*/
bool wasAdapted() const
{
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
-
- return (sumMarked != 0 || sumCoarsened != 0);
+ return (sumRefine_ != 0 || sumCoarsen_ != 0);
}
/*!
@@ -222,9 +219,6 @@ private:
template<class Indicator>
void adaptGrid(Indicator& indicator)
{
- // reset internal counter for marked elements
- marked_ = coarsened_ = 0;
-
// check for adaption interval: Adapt only at certain time step indices
if (problem_.timeManager().timeStepIndex() % adaptionInterval_ != 0)
return;
@@ -235,16 +229,17 @@ private:
indicator.calculateIndicator();
/**** 2) mark elements according to indicator *********/
- markElements(indicator);
+ int refine, coarsen;
+ std::tie(refine, coarsen) = markElements(indicator);
// abort if nothing in grid is marked
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
- if (sumMarked == 0 && sumCoarsened == 0)
+ sumRefine_ = problem_.grid().comm().sum(refine);
+ sumCoarsen_ = problem_.grid().comm().sum(coarsen);
+ if (sumRefine_ == 0 && sumCoarsen_ == 0)
return;
- else
- Dune::dinfo << marked_ << " cells have been marked_ to be refined, "
- << coarsened_ << " to be coarsened." << std::endl;
+ else if (problem_.grid().comm().rank() == 0)
+ Dune::dinfo << sumRefine_ << " cells have been marked to be refined, "
+ << sumCoarsen_ << " to be coarsened." << std::endl;
/**** 2b) Do pre-adaption step *****/
problem_.grid().preAdapt();
@@ -280,11 +275,10 @@ private:
* @param indicator The refinement indicator that is applied
*/
template<class Indicator>
- void markElements(Indicator& indicator)
+ std::pair<int, int> markElements(Indicator& indicator)
{
- typedef std::unordered_map<int, int> CoarsenMarkerType;
- CoarsenMarkerType coarsenMarker;
-
+ int refine = 0;
+ int coarsen = 0;
for (const auto& element : elements(problem_.gridView()))
{
// only mark non-ghost elements
@@ -294,18 +288,19 @@ private:
if (indicator.refine(element) && element.level() < levelMax_)
{
problem_.grid().mark( 1, element);
- ++marked_;
+ ++refine;
// this also refines the neighbor elements
- checkNeighborsRefine_(element);
+ refine += checkNeighborsRefine_(element);
}
if (indicator.coarsen(element) && element.hasFather())
{
problem_.grid().mark( -1, element );
- ++coarsened_;
+ ++coarsen;
}
}
}
+ return std::make_pair(refine, coarsen);
}
/*!
@@ -320,9 +315,10 @@ private:
* @param level level of the refined element: it is at least 1
* @return true if everything was successful
*/
- bool checkNeighborsRefine_(const Element &element, int level = 1)
+ int checkNeighborsRefine_(const Element &element, int level = 1)
{
// this also refines the neighbor elements
+ int refine = 0;
for(const auto& intersection : intersections(problem_.gridView(), element))
{
if(!intersection.neighbor())
@@ -338,13 +334,13 @@ private:
&& (outside.level() < element.level()))
{
problem_.grid().mark(1, outside);
- ++marked_;
+ ++refine;
if(level != levelMax_)
- checkNeighborsRefine_(outside, ++level);
+ refine += checkNeighborsRefine_(outside, ++level);
}
}
- return true;
+ return refine;
}
@@ -359,7 +355,7 @@ private:
*/
void forceRefineRatio(int maxLevelDelta = 1)
{
- LeafGridView leafGridView = problem_.gridView();
+ auto leafGridView = problem_.gridView();
// delete all existing marks
problem_.grid().postAdapt();
bool done;
@@ -394,7 +390,7 @@ private:
problem_.grid().postAdapt();
}
}
- while (done!=true);
+ while (!done);
}
// private Variables
@@ -402,8 +398,8 @@ private:
AdaptionHelper adaptionHelper_;
AdaptionIndicator adaptionIndicator_;
- int marked_;
- int coarsened_;
+ int sumRefine_;
+ int sumCoarsen_;
int levelMin_;
int levelMax_;
@@ -444,4 +440,4 @@ public:
};
}
-#endif /* DUMUX_IMPLICIT_GRIDADAPT_HH */
\ No newline at end of file
+#endif /* DUMUX_IMPLICIT_GRIDADAPT_HH */
diff --git a/dumux/implicit/adaptive/gridadaptinitializationindicator.hh b/dumux/implicit/adaptive/gridadaptinitializationindicator.hh
index b106447be6b16703bc81d256bcdbc208fb543e24..55c3db2bf2a0e90be0b2e3fe8177b5fe09a4a9b0 100644
--- a/dumux/implicit/adaptive/gridadaptinitializationindicator.hh
+++ b/dumux/implicit/adaptive/gridadaptinitializationindicator.hh
@@ -20,7 +20,6 @@
#define DUMUX_IMPLICIT_GRIDADAPTINITIALIZATIONINDICATOR_HH
#include <dune/geometry/type.hh>
-#include <dune/common/dynvector.hh>
#include "gridadaptproperties.hh"
/**
@@ -53,9 +52,9 @@ template<class TypeTag>
class ImplicitGridAdaptInitializationIndicator
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum
{
@@ -64,28 +63,23 @@ private:
numEq = GET_PROP_VALUE(TypeTag, NumEq),
};
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::Traits::template Codim<dim>::Entity Vertex;
- typedef typename GridView::Intersection Intersection;
- typedef typename GridView::Grid::ctype CoordScalar;
- typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
- typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
-
- typedef typename GET_PROP_TYPE(TypeTag, AdaptionIndicator) AdaptionIndicator;
-
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ using Element = typename GridView::Traits::template Codim<0>::Entity;
+ using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
+ using CoordScalar = typename GridView::Grid::ctype;
+ using ReferenceElement = typename Dune::ReferenceElement<CoordScalar, dim>;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+ using AdaptionIndicator = typename GET_PROP_TYPE(TypeTag, AdaptionIndicator);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-
enum { refineCell = 1 };
-
-
-
/*! \brief Search for a source term
*
* For every element we check if the element center or the element corners
@@ -97,10 +91,9 @@ private:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars)
{
- PrimaryVariables source(0.0);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; scvIdx++)
+ for (auto&& scv : scvs(fvGeometry))
{
- problem_.solDependentSource(source, element, fvGeometry, scvIdx, elemVolVars);
+ auto source = problem_.source(element, fvGeometry, elemVolVars, scv);
if (source.infinity_norm() > eps_)
return true;
}
@@ -118,57 +111,53 @@ private:
* \param element A grid element
* \param intersection The boundary intersection
*/
- bool hasRefineBC_(ElementBoundaryTypes &elemBcTypes,
- const Element& element,
- const Intersection& intersection,
+ bool hasRefineBC_(const Element& element,
const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars)
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementBoundaryTypes &elemBcTypes)
{
// Check for boundary conditions
if(isBox)
{
- Dune::GeometryType geoType = element.geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geoType);
- int fIdx = intersection.indexInInside();
- int numFaceVerts = refElement.size(fIdx, 1, dim);
- for (int faceVertexIdx = 0; faceVertexIdx < numFaceVerts; ++faceVertexIdx)
- {
- int scvIdx = refElement.subEntity(fIdx, 1, faceVertexIdx, dim);
- BoundaryTypes bcTypes = elemBcTypes[scvIdx];
- problemBoundaryTypes_(bcTypes, element.template subEntity<dim>(scvIdx));
- int bfIdx = fvGeometry.boundaryFaceIndex(fIdx, faceVertexIdx);
- for (int i = 0; i < numEq; i++)
- {
- if(bcTypes.isDirichlet(i) && refineAtDirichletBC_)
- return true;
- if(bcTypes.isNeumann(i) && refineAtFluxBC_)
- {
- PrimaryVariables fluxes(0.0);
- problem_.solDependentNeumann(fluxes, element, fvGeometry,
- intersection, scvIdx, bfIdx,
- elemVolVars);
- if (fluxes.infinity_norm() > eps_)
- return true;
- }
- }
- }
-
+ // Dune::GeometryType geoType = element.geometry().type();
+ // const ReferenceElement &refElement = ReferenceElements::general(geoType);
+ // int fIdx = intersection.indexInInside();
+ // int numFaceVerts = refElement.size(fIdx, 1, dim);
+ // for (int faceVertexIdx = 0; faceVertexIdx < numFaceVerts; ++faceVertexIdx)
+ // {
+ // int scvIdx = refElement.subEntity(fIdx, 1, faceVertexIdx, dim);
+ // BoundaryTypes bcTypes = elemBcTypes[scvIdx];
+ // problemBoundaryTypes_(bcTypes, element.template subEntity<dim>(scvIdx));
+ // int bfIdx = fvGeometry.boundaryFaceIndex(fIdx, faceVertexIdx);
+ // for (int i = 0; i < numEq; i++)
+ // {
+ // if(bcTypes.isDirichlet(i) && refineAtDirichletBC_)
+ // return true;
+ // if(bcTypes.isNeumann(i) && refineAtFluxBC_)
+ // {
+ // PrimaryVariables fluxes(0.0);
+ // problem_.solDependentNeumann(fluxes, element, fvGeometry,
+ // intersection, scvIdx, bfIdx,
+ // elemVolVars);
+ // if (fluxes.infinity_norm() > eps_)
+ // return true;
+ // }
+ // }
+ // }
+ return false;
}
else
{
- BoundaryTypes bcTypes = elemBcTypes[0];
- problem_.boundaryTypes(bcTypes, intersection);
- int bfIdx = intersection.indexInInside();
+ auto bcTypes = problem_.boundaryTypes(element, scvf);
for (int i = 0; i < numEq; i++)
{
if(bcTypes.isDirichlet(i) && refineAtDirichletBC_)
return true;
- if(bcTypes.isNeumann(i) && refineAtFluxBC_)
+
+ else if(bcTypes.isNeumann(i) && refineAtFluxBC_)
{
- PrimaryVariables fluxes(0.0);
- problem_.solDependentNeumann(fluxes, element, fvGeometry,
- intersection, 0, bfIdx,
- elemVolVars);
+ auto fluxes = problem_.neumann(element, fvGeometry, elemVolVars, scvf);
if (fluxes.infinity_norm() > eps_)
return true;
}
@@ -179,16 +168,19 @@ private:
// only actually call the problem method when it exists
template <class T = TypeTag>
- void problemBoundaryTypes_(BoundaryTypes& bcTypes,
- const typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox), Vertex>::type& v) const
+ BoundaryTypes problemBoundaryTypes_(const Element &element,
+ const typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox),
+ SubControlVolume>::type& scv) const
{
- problem_.boundaryTypes(bcTypes, v);
+ return problem_.boundaryTypes(element, scv);
}
template <class T = TypeTag>
- void problemBoundaryTypes_(BoundaryTypes& bcTypes,
- const typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox), Vertex>::type& v) const
- {}
-
+ BoundaryTypes problemBoundaryTypes_(const Element &element,
+ const typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox),
+ SubControlVolume>::type& scv) const
+ {
+ return BoundaryTypes();
+ }
public:
/*! \brief Calculates the indicator used for refinement/coarsening for each grid cell.
@@ -199,19 +191,17 @@ public:
if (!enableInitializationIndicator_)
return;
- //First adapt for boundary conditions and sources to get a good initial solution
+ // first adapt for boundary conditions and sources to get a good initial solution
if (nextMaxLevel_ == maxAllowedLevel_)
adaptionIndicator_.calculateIndicator();
// prepare an indicator for refinement
- indicatorVector_.resize(problem_.gridView().size(0));
- indicatorVector_ = 0;
+ indicatorVector_.resize(problem_.gridView().size(0), 0);
// 1) calculate Indicator -> min, maxvalues
// loop over all leaf elements
for (const auto& element : elements(problem_.gridView()))
{
-
int globalIdxI = problem_.elementMapper().index(element);
int level = element.level();
using std::max;
@@ -229,10 +219,11 @@ public:
continue;
// get the fvGeometry and elementVolVars needed for the bc and source interfaces
- FVElementGeometry fvGeometry;
- ElementVolumeVariables elemVolVars;
- fvGeometry.update(problem_.gridView(), element);
- elemVolVars.update(problem_, element, fvGeometry, false /* oldSol? */);
+ auto fvGeometry = localView(problem_.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(problem_.model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, problem_.model().curSol());
// Check if we have to refine around a source term
if (indicatorVector_[globalIdxI] != refineCell && refineAtSource_)
@@ -246,18 +237,18 @@ public:
}
// get the element boundary types
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_, element, fvGeometry);
+ ElementBoundaryTypes elemBcTypes;
+ elemBcTypes.update(problem_, element, fvGeometry);
// Check if we have to refine at the boundary
if (indicatorVector_[globalIdxI] != refineCell && (refineAtDirichletBC_ || refineAtFluxBC_))
{
// Calculate the boundary indicator for all boundary intersections
- for (const auto& intersection : intersections(problem_.gridView(), element))
+ for (auto&& scvf : scvfs(fvGeometry))
{
- if (intersection.boundary())
+ if (scvf.boundary())
{
- if(hasRefineBC_(bcTypes, element, intersection, fvGeometry, elemVolVars))
+ if(hasRefineBC_(element, fvGeometry, elemVolVars, scvf, elemBcTypes))
{
nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
indicatorVector_[globalIdxI] = refineCell;
@@ -275,7 +266,7 @@ public:
*
* \param element A grid element
*/
- bool refine(const Element& element)
+ bool refine(const Element& element) const
{
int idx = problem_.elementMapper().index(element);
@@ -293,28 +284,21 @@ public:
*
* \param element A grid element
*/
- bool coarsen(const Element& element)
- {
- return false;
- }
+ bool coarsen(const Element& element) const
+ { return false; }
- int maxLevel()
- {
- return maxLevel_;
- }
+ int maxLevel() const
+ { return maxLevel_; }
/*! \brief Initializes the adaption indicator class */
- void init()
- {};
+ void init() {}
/*! \brief If the model needs to be initialized after adaption.
* We always need to initialize since the hasRefineBC_ method needs information from
an up-to-date model.
*/
- bool initializeModel()
- {
- return true;
- }
+ constexpr bool initializeModel() const
+ { return true; }
/*! \brief Constructs a GridAdaptionIndicator instance
*
@@ -346,7 +330,7 @@ public:
private:
Problem& problem_;
AdaptionIndicator& adaptionIndicator_;
- Dune::DynamicVector<int> indicatorVector_;
+ std::vector<int> indicatorVector_;
int maxLevel_;
int nextMaxLevel_;
int minAllowedLevel_;
@@ -388,7 +372,7 @@ public:
*
* \param element A grid element
*/
- bool refine(const Element& element)
+ bool refine(const Element& element) const
{
return false;
}
@@ -399,19 +383,18 @@ public:
*
* \param element A grid element
*/
- bool coarsen(const Element& element)
+ bool coarsen(const Element& element) const
{
return false;
}
- bool initializeModel()
+ constexpr bool initializeModel() const
{
return false;
}
/*! \brief Initializes the adaption indicator class*/
- void init()
- {};
+ void init() {}
/*! \brief Constructs a GridAdaptionIndicator for initialization of an adaptive grid
*
diff --git a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
index 9f21d22dcd4673d4fe9f29cacd4d3333f9bb4833..80dc9e334f37cbc23f8e21f81728cb77179518e8 100644
--- a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
+++ b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
@@ -31,15 +31,17 @@ 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, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using ParentType = ImplicitAdaptionHelper<TypeTag>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum {
// Grid and world dimension
@@ -59,29 +61,25 @@ private:
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;
+ using Grid = typename GridView::Grid;
+ using Element = typename GridView::Traits::template Codim<0>::Entity;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename GridView::ctype CoordScalar;
- typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using CoordScalar = typename GridView::ctype;
+ using LocalPosition = Dune::FieldVector<CoordScalar, dim>;
- typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
- typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
+ using LocalFiniteElementCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using LocalFiniteElement = typename LocalFiniteElementCache::FiniteElementType;
struct AdaptedValues
{
- std::vector<PrimaryVariables> u;
+ ElementSolutionVector u;
int count;
PrimaryVariables associatedMass;
- AdaptedValues(): associatedMass(0)
- {
- count = 0;
- }
+ AdaptedValues(): count(0), associatedMass(0.0) {}
};
- typedef Dune::PersistentContainer<Grid, AdaptedValues> PersistentContainer;
+ using PersistentContainer = Dune::PersistentContainer<Grid, AdaptedValues>;
PersistentContainer adaptionMap_;
public:
@@ -92,7 +90,7 @@ public:
TwoPAdaptionHelper(Problem& problem) : ParentType(problem), adaptionMap_(problem.grid(), 0)
{
if(FluidSystem::isCompressible(wPhaseIdx) || FluidSystem::isCompressible(nPhaseIdx))
- DUNE_THROW(Dune::InvalidStateException, "Adaptionhelper is only for incompressible fluids mass-conservative!");
+ DUNE_THROW(Dune::InvalidStateException, "This adaption helper is only mass conservative for incompressible fluids!");
}
/*!
@@ -113,53 +111,43 @@ public:
// 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);
+ // get grid view on level grid
+ auto levelView = problem.grid().levelGridView(level);
for (const auto& element : elements(levelView))
{
//get your map entry
- AdaptedValues &adaptedValues = adaptionMap_[element];
+ auto& adaptedValues = adaptionMap_[element];
// put values in the map
if (element.isLeaf())
{
- FVElementGeometry fvGeometry;
- fvGeometry.update(problem.gridView(), element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- // get index
- int dofIdx = this->dofIndex(problem, element, scvIdx);
+ auto elemVolVars = localView(problem.model().curGlobalVolVars());
+ elemVolVars.bindElement(element, fvGeometry, problem.model().curSol());
- adaptedValues.u.push_back(problem.model().curSol()[dofIdx]);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ adaptedValues.u = problem.model().elementSolution(element, problem.model().curSol());
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);
+ volVars.update(adaptedValues.u, problem, element, scv);
+
+ adaptedValues.associatedMass[nPhaseIdx] += scv.volume() * volVars.density(nPhaseIdx)
+ * volVars.porosity() * volVars.saturation(nPhaseIdx);
+ adaptedValues.associatedMass[wPhaseIdx] += scv.volume() * volVars.density(wPhaseIdx)
+ * volVars.porosity() * volVars.saturation(wPhaseIdx);
}
adaptedValues.count = 1;
}
- //Average in father
+ // 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()];
- //For some grids the father element is identical to the son element.
- //For that case averaging is not necessary.
+ auto& adaptedValuesFather = adaptionMap_[element.father()];
+ // 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;
@@ -169,17 +157,8 @@ public:
if(isBox && !element.isLeaf())
{
- FVElementGeometry fvGeometry;
- fvGeometry.update(problem.gridView(), element);
-
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int dofIdx = this->dofIndex(problem, element, scvIdx);
-
- adaptedValues.u.push_back(problem.model().curSol()[dofIdx]);
- }
+ adaptedValues.u = problem.model().elementSolution(element, problem.model().curSol());
}
-
}
}
}
@@ -200,19 +179,19 @@ public:
void reconstructPrimVars(Problem& problem)
{
adaptionMap_.resize();
- //vectors storing the mass associated with each vertex, when using the box method
+ // vectors storing the mass associated with each vertex, when using the box method
std::vector<Scalar> massCoeff;
std::vector<Scalar> associatedMass;
if(isBox)
{
- massCoeff.resize(problem.model().numDofs(),0.0);
- associatedMass.resize(problem.model().numDofs(),0.0);
+ 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);
+ auto levelView = problem.grid().levelGridView(level);
for (const auto& element : elements(levelView))
{
@@ -222,54 +201,45 @@ public:
if (!element.isNew() || element.level() == 0)
{
- //entry is in map, write in leaf
+ // entry is in map, write in leaf
if (element.isLeaf())
{
- AdaptedValues &adaptedValues = adaptionMap_[element];
+ auto& adaptedValues = adaptionMap_[element];
- FVElementGeometry fvGeometry;
- fvGeometry.update(problem.gridView(), element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- // get index
- int dofIdx = this->dofIndex(problem, element, scvIdx);
+ auto elementVolume = element.geometry().volume();
+ for (auto&& scv : scvs(fvGeometry))
+ {
VolumeVariables volVars;
- volVars.update(adaptedValues.u[scvIdx],
- problem,
- element,
- fvGeometry,
- scvIdx,
- false);
-
- Scalar volumeElement = fvGeometry.elementVolume;
+ volVars.update(adaptedValues.u, problem, element, scv);
- this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdx],scvIdx);
+ problem.model().curSol()[scv.dofIndex()] = getPriVars(adaptedValues, scv);
- if (int(formulation) == pwsn)
+ if (formulation == pwsn)
{
- problem.model().curSol()[dofIdx][saturationIdx] = adaptedValues.associatedMass[nPhaseIdx];
- problem.model().curSol()[dofIdx][saturationIdx] /= volumeElement * volVars.density(nPhaseIdx) * volVars.porosity();
+ problem.model().curSol()[scv.dofIndex()][saturationIdx] = adaptedValues.associatedMass[nPhaseIdx];
+ problem.model().curSol()[scv.dofIndex()][saturationIdx] /= elementVolume * volVars.density(nPhaseIdx) * volVars.porosity();
}
- else if (int(formulation) == pnsw)
+ else if (formulation == pnsw)
{
- problem.model().curSol()[dofIdx][saturationIdx] = adaptedValues.associatedMass[wPhaseIdx];
- problem.model().curSol()[dofIdx][saturationIdx] /= volumeElement * volVars.density(wPhaseIdx) * volVars.porosity();
+ problem.model().curSol()[scv.dofIndex()][saturationIdx] = adaptedValues.associatedMass[wPhaseIdx];
+ problem.model().curSol()[scv.dofIndex()][saturationIdx] /= elementVolume * volVars.density(wPhaseIdx) * volVars.porosity();
}
if(isBox)
{
- Scalar volume = fvGeometry.subContVol[scvIdx].volume;
if (int(formulation) == pwsn)
{
- massCoeff[dofIdx] += volume * volVars.density(nPhaseIdx) * volVars.porosity();
- associatedMass[dofIdx] += volume/volumeElement*adaptedValues.associatedMass[nPhaseIdx];
+ massCoeff[scv.dofIndex()] += scv.volume() * volVars.density(nPhaseIdx) * volVars.porosity();
+ associatedMass[scv.dofIndex()] += scv.volume() / elementVolume * adaptedValues.associatedMass[nPhaseIdx];
}
else if (int(formulation) == pnsw)
{
- massCoeff[dofIdx] += volume * volVars.density(wPhaseIdx) * volVars.porosity();
- associatedMass[dofIdx] += volume/volumeElement*adaptedValues.associatedMass[wPhaseIdx];
+ massCoeff[scv.dofIndex()] += scv.volume() * volVars.density(wPhaseIdx) * volVars.porosity();
+ associatedMass[scv.dofIndex()] += scv.volume() / elementVolume * adaptedValues.associatedMass[wPhaseIdx];
}
}
@@ -282,127 +252,118 @@ public:
// value is not in map, interpolate from father element
if (element.hasFather())
{
- auto eFather = element.father();
- while(eFather.isNew() && eFather.level() > 0)
- eFather = eFather.father();
+ auto fatherElement = element.father();
+ while(fatherElement.isNew() && fatherElement.level() > 0)
+ fatherElement = fatherElement.father();
Scalar massFather = 0.0;
if(!isBox)
{
- AdaptedValues& adaptedValuesFather = adaptionMap_[eFather];
+ auto& adaptedValuesFather = adaptionMap_[fatherElement];
- if (int(formulation) == pwsn)
- {
+ if (formulation == pwsn)
massFather = adaptedValuesFather.associatedMass[nPhaseIdx];
- }
- else if (int(formulation) == pnsw)
- {
+
+ else if (formulation == pnsw)
massFather = adaptedValuesFather.associatedMass[wPhaseIdx];
- }
// access new son
- AdaptedValues& adaptedValues = adaptionMap_[element];
+ auto& adaptedValues = adaptionMap_[element];
adaptedValues.count = 1;
- FVElementGeometry fvGeometry;
- fvGeometry.update(problem.gridView(), element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- adaptedValues.u.push_back(adaptedValuesFather.u[0]);
+ adaptedValues.u = adaptedValuesFather.u;
- 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)
+ for (auto&& scv : scvs(fvGeometry))
{
- massCoeffSon = volume * volVars.density(wPhaseIdx) * volVars.porosity();
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u, problem, element, scv);
+
+ Scalar massCoeffSon = 0.0;
+ if (int(formulation) == pwsn)
+ massCoeffSon = scv.volume() * volVars.density(nPhaseIdx) * volVars.porosity();
+
+ else if (int(formulation) == pnsw)
+ massCoeffSon = scv.volume() * volVars.density(wPhaseIdx) * volVars.porosity();
+
+ auto fatherElementVolume = fatherElement.geometry().volume();
+ adaptedValues.u[0][saturationIdx] = (scv.volume() / fatherElementVolume * massFather)/massCoeffSon;
}
- 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())
{
// access new CellData object
- int newIdxI = this->elementIndex(problem, element);
-
- this->setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI],0);
+ auto newIdxI = problem.elementMapper().index(element);
+ problem.model().curSol()[newIdxI] = adaptedValues.u[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();
-
- // Interpolate values from father element by using ansatz functions
- 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, element, scvIdx);
- 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);
- }
-
- }
+ // auto& adaptedValuesFather = adaptionMap_[fatherElement];
+ // // access new son
+ // auto& 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 = fatherElement.geometry().type();
+ //
+ // // Interpolate values from father element by using ansatz functions
+ // 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 fatherElementVolume = fatherElement.geometry().volume();
+ //
+ // int dofIdxGlobal = this->dofIndex(problem, element, scvIdx);
+ // if (int(formulation) == pwsn)
+ // {
+ // massCoeff[dofIdxGlobal] += volume * volVars.density(nPhaseIdx) * volVars.porosity();
+ // associatedMass[dofIdxGlobal] += volume/fatherElementVolume*adaptedValuesFather.associatedMass[nPhaseIdx];
+ // }
+ // else if (int(formulation) == pnsw)
+ // {
+ // massCoeff[dofIdxGlobal] += volume * volVars.density(wPhaseIdx) * volVars.porosity();
+ // associatedMass[dofIdxGlobal] += volume/fatherElementVolume*adaptedValuesFather.associatedMass[wPhaseIdx];
+ // }
+ //
+ // this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdxGlobal],scvIdx);
+ // }
+ //
+ // }
}
}
else
@@ -417,10 +378,8 @@ public:
if(isBox)
{
- for(int dofIdx = 0; dofIdx < problem.model().numDofs(); dofIdx++)
- {
- problem.model().curSol()[dofIdx][saturationIdx] = associatedMass[dofIdx]/massCoeff[dofIdx];
- }
+ for(int dofIdxGlobal = 0; dofIdxGlobal < problem.model().numDofs(); dofIdxGlobal++)
+ problem.model().curSol()[dofIdxGlobal][saturationIdx] = associatedMass[dofIdxGlobal] / massCoeff[dofIdxGlobal];
}
// reset entries in restrictionmap
@@ -454,9 +413,6 @@ public:
{
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;
@@ -475,21 +431,21 @@ public:
* \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)
+ static PrimaryVariables getPriVars(const AdaptedValues& adaptedValues, const SubControlVolume& scv)
{
- PrimaryVariables uNew(0);
if(!isBox)
{
- uNew = adaptedValues.u[scvIdx];
+ auto uNew = adaptedValues.u[0];
uNew /= adaptedValues.count;
+ return uNew;
}
else
{
- uNew = adaptedValues.u[scvIdx];
+ return adaptedValues.u[scv.index()];
}
-
- u = uNew;
}
};
-}
+
+} // end namespace Dumux
+
#endif