From aee63233791fa0e9dfc7eb16d29add59ca1bae33 Mon Sep 17 00:00:00 2001 From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de> Date: Tue, 21 Apr 2015 07:20:54 +0000 Subject: [PATCH] Implementation of grid adaptation for box scheme with conforming grid refinement. Renamed classes in a consistent way. Reviewed by Timo git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/gridadapt@14600 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- ...indicator2p.hh => gridadaptindicator2p.hh} | 92 ++++++- dumux/implicit/Makefile.am | 3 +- dumux/implicit/adaptive/Makefile.am | 4 + .../{common => adaptive}/adaptationhelper.hh | 237 +++++++++++++----- .../{common => adaptive}/gridadapt.hh | 50 ++-- .../gridadaptindicatordefault.hh | 10 +- .../gridadaptinitializationindicator.hh | 38 +-- .../gridadaptproperties.hh | 12 +- .../gridadaptpropertydefaults.hh | 10 +- dumux/implicit/common/implicitmodel.hh | 2 +- dumux/implicit/common/implicitproblem.hh | 6 +- dumux/implicit/common/implicitproperties.hh | 2 +- 12 files changed, 334 insertions(+), 132 deletions(-) rename dumux/implicit/2p/{gridadaptionindicator2p.hh => gridadaptindicator2p.hh} (68%) create mode 100644 dumux/implicit/adaptive/Makefile.am rename dumux/implicit/{common => adaptive}/adaptationhelper.hh (53%) rename dumux/implicit/{common => adaptive}/gridadapt.hh (91%) rename dumux/implicit/{common => adaptive}/gridadaptindicatordefault.hh (88%) rename dumux/implicit/{common => adaptive}/gridadaptinitializationindicator.hh (91%) rename dumux/implicit/{common => adaptive}/gridadaptproperties.hh (91%) rename dumux/implicit/{common => adaptive}/gridadaptpropertydefaults.hh (87%) diff --git a/dumux/implicit/2p/gridadaptionindicator2p.hh b/dumux/implicit/2p/gridadaptindicator2p.hh similarity index 68% rename from dumux/implicit/2p/gridadaptionindicator2p.hh rename to dumux/implicit/2p/gridadaptindicator2p.hh index a2b798fcfa..35b9a260f7 100644 --- a/dumux/implicit/2p/gridadaptionindicator2p.hh +++ b/dumux/implicit/2p/gridadaptindicator2p.hh @@ -16,27 +16,28 @@ * 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_IMPLICIT_GRIDADAPTIONINDICATOR2P_HH -#define DUMUX_IMPLICIT_GRIDADAPTIONINDICATOR2P_HH +#ifndef DUMUX_IMPLICIT_GRIDADAPTINDICATOR2P_HH +#define DUMUX_IMPLICIT_GRIDADAPTINDICATOR2P_HH #include <dune/common/version.hh> #include "2pproperties.hh" +#include <dune/localfunctions/lagrange/pqkfactory.hh> //#include <dumux/linear/vectorexchange.hh> /** * @file - * @brief Class defining a standard, saturation dependent indicator for grid adaption + * @brief Class defining a standard, saturation dependent indicator for grid adaptation */ namespace Dumux { /*!\ingroup IMPES - * @brief Class defining a standard, saturation dependent indicator for grid adaption + * @brief Class defining a standard, saturation dependent indicator for grid adaptation * * \tparam TypeTag The problem TypeTag */ template<class TypeTag> -class ImplicitGridAdaptionIndicator2P +class ImplicitGridAdaptIndicator2P { private: typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; @@ -60,6 +61,21 @@ private: nPhaseIdx = Indices::nPhaseIdx }; + enum { + // Grid and world dimension + dim = GridView::dimension, + dimWorld = GridView::dimensionworld + }; + + enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + enum { dofCodim = isBox ? dim : 0 }; + + typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + typedef typename GridView::ctype CoordScalar; + + typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache; + typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement; + public: /*! \brief Calculates the indicator used for refinement/coarsening for each grid cell. * @@ -68,9 +84,9 @@ public: void calculateIndicator() { // prepare an indicator for refinement - if(indicatorVector_.size() != problem_.model().numDofs()) + if(indicatorVector_.size() != problem_.gridView().size(0)) { - indicatorVector_.resize(problem_.model().numDofs()); + indicatorVector_.resize(problem_.gridView().size(0)); } indicatorVector_ = -1e100; @@ -91,7 +107,31 @@ public: int globalIdxI = problem_.elementMapper().map(*eIt); #endif - Scalar satI = problem_.model().curSol()[globalIdxI][saturationIdx]; + Scalar satI = 0.0; + + if(!isBox) + satI = problem_.model().curSol()[globalIdxI][saturationIdx]; + else + { + const LocalFiniteElementCache feCache; + const auto geometryI = eIt->geometry(); + Dune::GeometryType geomType = geometryI.type(); + + GlobalPosition centerI = geometryI.local(geometryI.center()); + const LocalFiniteElement &localFiniteElement = feCache.get(geomType); + std::vector<Dune::FieldVector<Scalar, 1> > shapeVal; + localFiniteElement.localBasis().evaluateFunction(centerI, shapeVal); + + for (int i = 0; i < shapeVal.size(); ++i) + { +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) + int dofIdxGlobal = problem_.model().dofMapper().subIndex(*eIt, i, dofCodim); +#else + int dofIdxGlobal = problem_.model().dofMapper().map(*eIt, i, dofCodim); +#endif + satI += shapeVal[i]*problem_.model().curSol()[dofIdxGlobal][saturationIdx]; + } + } globalMin = std::min(satI, globalMin); globalMax = std::max(satI, globalMax); @@ -115,7 +155,33 @@ public: // Visit intersection only once if (eIt->level() > outside->level() || (eIt->level() == outside->level() && globalIdxI < globalIdxJ)) { - Scalar satJ = problem_.model().curSol()[globalIdxJ][saturationIdx]; + Scalar satJ = 0.0; + + if(!isBox) + satJ = problem_.model().curSol()[globalIdxJ][saturationIdx]; + else + { + const LocalFiniteElementCache feCache; + const auto geometryJ = outside->geometry(); + Dune::GeometryType geomType = geometryJ.type(); + + GlobalPosition centerJ = geometryJ.local(geometryJ.center()); + const LocalFiniteElement &localFiniteElement = feCache.get(geomType); + std::vector<Dune::FieldVector<Scalar, 1> > shapeVal; + localFiniteElement.localBasis().evaluateFunction(centerJ, shapeVal); + + for (int i = 0; i < shapeVal.size(); ++i) + { +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) + int dofIdxGlobal = problem_.model().dofMapper().subIndex(*outside, i, dofCodim); +#else + int dofIdxGlobal = problem_.model().dofMapper().map(*outside, i, dofCodim); +#endif + satJ += shapeVal[i]*problem_.model().curSol()[dofIdxGlobal][saturationIdx]; + } + } + + Scalar localdelta = std::abs(satI - satJ); indicatorVector_[globalIdxI][0] = std::max(indicatorVector_[globalIdxI][0], localdelta); @@ -174,15 +240,15 @@ public: #endif } - /*! \brief Initializes the adaption indicator class*/ + /*! \brief Initializes the adaptation indicator class*/ void init() { refineBound_ = 0.; coarsenBound_ = 0.; - indicatorVector_.resize(problem_.model().numDofs()); + indicatorVector_.resize(problem_.gridView().size(0)); }; - /*! @brief Constructs a GridAdaptionIndicator instance + /*! @brief Constructs a GridAdaptIndicator instance * * This standard indicator is based on the saturation gradient. * It checks the local gradient compared to the maximum global gradient. @@ -191,7 +257,7 @@ public: * * \param problem The problem object */ - ImplicitGridAdaptionIndicator2P (Problem& problem): + ImplicitGridAdaptIndicator2P (Problem& problem): problem_(problem) { refinetol_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, GridAdapt, RefineTolerance); diff --git a/dumux/implicit/Makefile.am b/dumux/implicit/Makefile.am index 68d6a4c1f3..b2c84a0204 100644 --- a/dumux/implicit/Makefile.am +++ b/dumux/implicit/Makefile.am @@ -11,7 +11,8 @@ SUBDIRS = 1p \ common \ mpnc \ nonisothermal \ - richards + richards \ + adaptive implicitdir = $(includedir)/dumux/implicit diff --git a/dumux/implicit/adaptive/Makefile.am b/dumux/implicit/adaptive/Makefile.am new file mode 100644 index 0000000000..9051d7c418 --- /dev/null +++ b/dumux/implicit/adaptive/Makefile.am @@ -0,0 +1,4 @@ +adaptivedir = $(includedir)/dumux/implicit/adaptive +adaptive_HEADERS := $(wildcard *.hh) + +include $(top_srcdir)/am/global-rules diff --git a/dumux/implicit/common/adaptationhelper.hh b/dumux/implicit/adaptive/adaptationhelper.hh similarity index 53% rename from dumux/implicit/common/adaptationhelper.hh rename to dumux/implicit/adaptive/adaptationhelper.hh index 8d692692bf..1a1d7539aa 100644 --- a/dumux/implicit/common/adaptationhelper.hh +++ b/dumux/implicit/adaptive/adaptationhelper.hh @@ -20,6 +20,7 @@ #define DUMUX_ADAPTATIONHELPER_HH #include <dune/grid/utility/persistentcontainer.hh> +#include <dune/localfunctions/lagrange/pqkfactory.hh> //#include <dumux/linear/vectorexchange.hh> /** @@ -37,6 +38,7 @@ 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; struct AdaptedValues { @@ -48,13 +50,32 @@ private: } }; + 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::Grid Grid; typedef typename Grid::LevelGridView LevelGridView; - typedef typename LevelGridView::template Codim<0>::Iterator LevelIterator; + typedef typename LevelGridView::template Codim<dofCodim>::Iterator LevelIterator; + typedef typename LevelGridView::template Codim<0>::Iterator ElementLevelIterator; typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer; typedef typename GridView::Traits::template Codim<0>::Entity Element; + typedef typename GridView::Traits::template Codim<dofCodim>::Entity DofEntity; + typedef typename GridView::Traits::template Codim<dofCodim>::EntityPointer DofPointer; typedef Dune::PersistentContainer<Grid, AdaptedValues> PersistentContainer; + 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; + private: const GridView gridView_; const Grid& grid_; @@ -69,7 +90,7 @@ public: * @param gridView a DUNE gridview object corresponding to diffusion and transport equation */ AdaptationHelper(const GridView& gridView) : - gridView_(gridView), grid_(gridView.grid()), adaptationMap_(grid_, 0) + gridView_(gridView), grid_(gridView.grid()), adaptationMap_(grid_, dofCodim) {} @@ -94,31 +115,51 @@ public: //get grid view on level grid LevelGridView levelView = grid_.levelGridView(level); - for (LevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt) + if(!isBox) + { + for (ElementLevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt) + { + //get your map entry + AdaptedValues &adaptedValues = adaptationMap_[*eIt]; + + // put your value in the map + if (eIt->isLeaf()) + { + // get index + int indexI = this->elementIndex(problem, *eIt); + + storeAdaptationValues(adaptedValues, problem.model().curSol()[indexI]); + + adaptedValues.count = 1; + } + //Average in father + if (eIt->level() > 0) + { + ElementPointer epFather = eIt->father(); + int indexI = this->elementIndex(problem, *epFather); + AdaptedValues& adaptedValuesFather = adaptationMap_[*epFather]; + adaptedValuesFather.count += 1; + storeAdaptationValues(adaptedValues, adaptedValuesFather, + problem.model().curSol()[indexI]); + } + + } + } + else { - //get your map entry - AdaptedValues &adaptedValues = adaptationMap_[*eIt]; + for (LevelIterator dofIt = levelView.template begin<dofCodim>(); dofIt != levelView.template end<dofCodim>(); ++dofIt) + { + //get your map entry + AdaptedValues &adaptedValues = adaptationMap_[*dofIt]; - // put your value in the map - if (eIt->isLeaf()) - { - // get index - int indexI = this->elementIndex(*eIt, problem); + // put your value in the map + int indexI = this->dofIndex(problem, *dofIt); - storeAdaptionValues(adaptedValues, problem.model().curSol()[indexI]); + storeAdaptationValues(adaptedValues, problem.model().curSol()[indexI]); - adaptedValues.count = 1; - } - //Average in father - if (eIt->level() > 0) - { - ElementPointer epFather = eIt->father(); - int indexI = this->elementIndex(*epFather, problem); - AdaptedValues& adaptedValuesFather = adaptationMap_[*epFather]; - adaptedValuesFather.count += 1; - storeAdaptionValues(adaptedValues, adaptedValuesFather, - problem.model().curSol()[indexI]); - } + adaptedValues.count = 1; + + } } } } @@ -144,7 +185,7 @@ public: { LevelGridView levelView = grid_.levelGridView(level); - for (LevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt) + for (ElementLevelIterator eIt = levelView.template begin<0>(); eIt != levelView.template end<0>(); ++eIt) { // only treat non-ghosts, ghost data is communicated afterwards if (eIt->partitionType() == Dune::GhostEntity) @@ -155,36 +196,103 @@ public: //entry is in map, write in leaf if (eIt->isLeaf()) { - AdaptedValues &adaptedValues = adaptationMap_[*eIt]; - int newIdxI = this->elementIndex(*eIt, problem); + if(!isBox) + { + AdaptedValues &adaptedValues = adaptationMap_[*eIt]; + int newIdxI = this->elementIndex(problem, *eIt); + + setAdaptationValues(adaptedValues, problem.model().curSol()[newIdxI]); + } + else + { +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) + unsigned int numSubEntities = eIt->subEntities(dofCodim); +#else + int numSubEntities = eIt->template count <dofCodim>(); +#endif + + for(unsigned int i = 0; i < numSubEntities; i++) + { + DofPointer subEntity = eIt->template subEntity <dofCodim>(i); + AdaptedValues &adaptedValues = adaptationMap_[*subEntity]; + int newIdxI = this->dofIndex(problem, *subEntity); - setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]); + setAdaptationValues(adaptedValues, problem.model().curSol()[newIdxI]); + + } + } } } else { // value is not in map, interpolate from father element - if (eIt->level() > 0) + if (eIt->level() > 0 && eIt->hasFather()) { ElementPointer epFather = eIt->father(); - // create new entry: reconstruct from adaptationMap_[*father] to a new - // adaptationMap_[*son] - reconstructAdaptionValues(adaptationMap_, *epFather, *eIt, problem); - - // access new son - AdaptedValues& adaptedValues = adaptationMap_[*eIt]; - adaptedValues.count = 1; - - // if we are on leaf, store reconstructed values of son in CellData object - if (eIt->isLeaf()) + if(!isBox) { - // acess new CellData object - int newIdxI = this->elementIndex(*eIt, problem); - - setAdaptionValues(adaptedValues, problem.model().curSol()[newIdxI]); + // create new entry: reconstruct from adaptationMap_[*father] to a new + // adaptationMap_[*son] + reconstructAdaptationValues(adaptationMap_, *epFather, *eIt, problem); + + // access new son + AdaptedValues& adaptedValues = adaptationMap_[*eIt]; + adaptedValues.count = 1; + + // if we are on leaf, store reconstructed values of son in CellData object + if (eIt->isLeaf()) + { + // acess new CellData object + int newIdxI = this->elementIndex(problem, *eIt); + + setAdaptationValues(adaptedValues, problem.model().curSol()[newIdxI]); + } + } + else + { +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) + unsigned int numSubEntities= eIt->subEntities(dofCodim); +#else + int numSubEntities= eIt->template count <dofCodim>(); +#endif + const auto geometryI = eIt->geometry(); + + for(unsigned int i = 0; i < numSubEntities; i++) + { + DofPointer subEntity = eIt->template subEntity <dofCodim>(i); + AdaptedValues &adaptedValues = adaptationMap_[*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) + { + DofPointer subEntityFather = epFather->template subEntity <dofCodim>(j); + AdaptedValues & adaptedValuesFather = adaptationMap_[*subEntityFather]; + u.axpy(shapeVal[j], adaptedValuesFather.u); + } + + adaptedValues.u = u; + adaptedValues.count = 1; + } + + if (eIt->isLeaf()) + { + int newIdxI = this->dofIndex(problem, *subEntity); + setAdaptationValues(adaptedValues, problem.model().curSol()[newIdxI]); + } + + } } } + } } @@ -214,7 +322,7 @@ public: * \param adaptedValues Container for model-specific values to be adapted * \param element The element to be stored */ - static void storeAdaptionValues(AdaptedValues& adaptedValues, const PrimaryVariables& u) + static void storeAdaptationValues(AdaptedValues& adaptedValues, const PrimaryVariables& u) { adaptedValues.u = u; } @@ -228,12 +336,19 @@ public: * \param adaptedValuesFather Values to be adapted of father cell * \param fatherElement The element of the father */ - static void storeAdaptionValues(AdaptedValues& adaptedValues, + static void storeAdaptationValues(AdaptedValues& adaptedValues, AdaptedValues& adaptedValuesFather, const PrimaryVariables& u) { - adaptedValuesFather.u += adaptedValues.u; - adaptedValuesFather.u /= adaptedValues.count; + if(!isBox) + { + adaptedValuesFather.u += adaptedValues.u; + adaptedValuesFather.u /= adaptedValues.count; + } + else + { + adaptedValuesFather.u = adaptedValues.u; + } } //! Set adapted values in CellData /** @@ -244,7 +359,7 @@ public: * \param adaptedValues Container for model-specific values to be adapted * \param element The element where things are stored. */ - static void setAdaptionValues(AdaptedValues& adaptedValues, PrimaryVariables& u) + static void setAdaptationValues(AdaptedValues& adaptedValues, PrimaryVariables& u) { PrimaryVariables uNew = adaptedValues.u; uNew /= adaptedValues.count; @@ -258,28 +373,36 @@ public: * generated son cell. New cell is stored into the global * adaptationMap. * - * \param adaptionMap Global map storing all values to be adapted + * \param adaptationMap 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) + static void reconstructAdaptationValues(Dune::PersistentContainer<Grid, AdaptedValues>& adaptationMap, + const Element& father, const Element& son, const Problem& problem) { - AdaptedValues& adaptedValues = adaptionMap[son]; - AdaptedValues& adaptedValuesFather = adaptionMap[father]; + AdaptedValues& adaptedValues = adaptationMap[son]; + AdaptedValues& adaptedValuesFather = adaptationMap[father]; - adaptedValues.u = adaptedValuesFather.u; - adaptedValues.u /= adaptedValuesFather.count; + adaptedValues.u = adaptedValuesFather.u; + adaptedValues.u /= adaptedValuesFather.count; + } + int dofIndex(const Problem& problem, const DofEntity& entity) const + { +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) + return problem.model().dofMapper().index(entity); +#else + return problem.model().dofMapper().map(entity); +#endif } - int elementIndex(const Element& element, const Problem& problem) const + int elementIndex(const Problem& problem, const Element& element) const { #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4) - return problem.elementMapper().index(element); + return problem.elementMapper().index(element); #else - return problem.elementMapper().map(element); + return problem.elementMapper().map(element); #endif } diff --git a/dumux/implicit/common/gridadapt.hh b/dumux/implicit/adaptive/gridadapt.hh similarity index 91% rename from dumux/implicit/common/gridadapt.hh rename to dumux/implicit/adaptive/gridadapt.hh index c154801726..4d124df513 100644 --- a/dumux/implicit/common/gridadapt.hh +++ b/dumux/implicit/adaptive/gridadapt.hh @@ -56,8 +56,8 @@ class ImplicitGridAdapt typedef typename Grid::template Codim<0>::Entity Element; typedef typename Grid::template Codim<0>::EntityPointer ElementPointer; - typedef typename GET_PROP_TYPE(TypeTag, AdaptionIndicator) AdaptionIndicator; - typedef typename GET_PROP_TYPE(TypeTag, AdaptionInitializationIndicator) AdaptionInitializationIndicator; + typedef typename GET_PROP_TYPE(TypeTag, AdaptationIndicator) AdaptationIndicator; + typedef typename GET_PROP_TYPE(TypeTag, AdaptationInitializationIndicator) AdaptationInitializationIndicator; public: /*! @@ -65,11 +65,11 @@ public: * @param problem The problem */ ImplicitGridAdapt (Problem& problem) - : problem_(problem), adaptationHelper_(problem.gridView()), adaptionIndicator_(problem), marked_(0), coarsened_(0) + : problem_(problem), adaptationHelper_(problem.gridView()), adaptationIndicator_(problem), marked_(0), coarsened_(0) { levelMin_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MinLevel); levelMax_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MaxLevel); - adaptationInterval_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, AdaptionInterval); + adaptationInterval_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, AdaptationInterval); if (levelMin_ < 0) { @@ -82,29 +82,29 @@ public: * * Prepares the grid for simulation after the initialization of the * problem. The applied indicator is selectable via the property - * AdaptionInitializationIndicator + * AdaptationInitializationIndicator */ void init() { - adaptionIndicator_.init(); + adaptationIndicator_.init(); if (!GET_PARAM_FROM_GROUP(TypeTag, bool, GridAdapt, EnableInitializationIndicator)) return; - AdaptionInitializationIndicator adaptionInitIndicator(problem_, adaptionIndicator_); + AdaptationInitializationIndicator adaptationInitIndicator(problem_, adaptationIndicator_); int maxIter = 2*levelMax_; int iter = 0; while (iter <= maxIter) { - adaptGrid(adaptionInitIndicator); + adaptGrid(adaptationInitIndicator); if (!wasAdapted()) { break; } - int shouldInitialize = adaptionInitIndicator.initializeModel(); + int shouldInitialize = adaptationInitIndicator.initializeModel(); if (problem_.grid().comm().max(shouldInitialize)) { problem_.model().init(problem_); @@ -119,7 +119,7 @@ public: * * This method is called from IMPETProblem::preTimeStep() if * adaptive grids are used in the simulation. It uses the standard - * indicator (selected by the property AdaptionIndicator) and forwards to + * indicator (selected by the property AdaptationIndicator) and forwards to * with it to the ultimate method adaptGrid(indicator), which * uses a standard procedure for adaptivity: * 1) Determine the refinement indicator @@ -130,7 +130,7 @@ public: */ void adaptGrid() { - adaptGrid(adaptionIndicator_) ; + adaptGrid(adaptationIndicator_) ; } /*! @@ -155,7 +155,7 @@ public: // reset internal counter for marked elements marked_ = coarsened_ = 0; - // check for adaption interval: Adapt only at certain time step indices + // check for adaptation interval: Adapt only at certain time step indices if (problem_.timeManager().timeStepIndex() % adaptationInterval_ != 0) return; @@ -176,7 +176,10 @@ public: Dune::dinfo << marked_ << " cells have been marked_ to be refined, " << coarsened_ << " to be coarsened." << std::endl; - /**** 2b) Do pre-adaption step *****/ + std::cout << marked_ << " cells have been marked_ to be refined, " + << coarsened_ << " to be coarsened." << std::endl; + + /**** 2b) Do pre-adaptation step *****/ problem_.grid().preAdapt(); problem_.preAdapt(); @@ -188,6 +191,7 @@ public: // update mapper to new cell indices problem_.elementMapper().update(); + problem_.vertexMapper().update(); // adapt size of vectors problem_.model().adaptVariableSize(); @@ -233,7 +237,9 @@ public: } if (indicator.coarsen(*eIt) && eIt->hasFather()) { - int idx = idSet.id(*(eIt->father())); + problem_.grid().mark( -1, *eIt ); + ++coarsened_; + /*int idx = idSet.id(*(eIt->father())); typename CoarsenMarkerType::iterator it = coarsenMarker.find(idx); if (it != coarsenMarker.end()) { @@ -242,11 +248,11 @@ public: else { coarsenMarker[idx] = 1; - } + }*/ } } // coarsen - for (LeafIterator eIt = problem_.gridView().template begin<0>(); + /*for (LeafIterator eIt = problem_.gridView().template begin<0>(); eIt!=problem_.gridView().template end<0>(); ++eIt) { // only mark non-ghost elements @@ -286,7 +292,7 @@ public: } } } - } + }*/ } /*! @@ -337,14 +343,14 @@ public: return levelMin_; } - AdaptionIndicator& adaptionIndicator() + AdaptationIndicator& adaptationIndicator() { - return adaptionIndicator_; + return adaptationIndicator_; } - AdaptionIndicator& adaptionIndicator() const + AdaptationIndicator& adaptationIndicator() const { - return adaptionIndicator_; + return adaptationIndicator_; } private: @@ -448,7 +454,7 @@ private: // private Variables Problem& problem_; - AdaptionIndicator adaptionIndicator_; + AdaptationIndicator adaptationIndicator_; int marked_; int coarsened_; diff --git a/dumux/implicit/common/gridadaptindicatordefault.hh b/dumux/implicit/adaptive/gridadaptindicatordefault.hh similarity index 88% rename from dumux/implicit/common/gridadaptindicatordefault.hh rename to dumux/implicit/adaptive/gridadaptindicatordefault.hh index 399d626825..331892e18f 100644 --- a/dumux/implicit/common/gridadaptindicatordefault.hh +++ b/dumux/implicit/adaptive/gridadaptindicatordefault.hh @@ -21,12 +21,12 @@ /** * @file - * @brief Class defining a default indicator for grid adaption + * @brief Class defining a default indicator for grid adaptation */ namespace Dumux { /*!\ingroup ImplicitGridAdaptIndicator - * @brief Class defining a default indicator for grid adaption + * @brief Class defining a default indicator for grid adaptation * *Default implementation * @@ -69,16 +69,16 @@ public: return false; } - /*! \brief Initializes the adaption indicator class*/ + /*! \brief Initializes the adaptation indicator class*/ void init() {}; - /*! \brief Constructs a GridAdaptionIndicator for initialization of an adaptive grid + /*! \brief Constructs a GridAdaptationIndicator for initialization of an adaptive grid * * Default implementation * * \param problem The problem object - * \param adaptionIndicator Indicator whether a be adapted + * \param adaptationIndicator Indicator whether a be adapted */ ImplicitGridAdaptIndicatorDefault(Problem& problem) {} diff --git a/dumux/implicit/common/gridadaptinitializationindicator.hh b/dumux/implicit/adaptive/gridadaptinitializationindicator.hh similarity index 91% rename from dumux/implicit/common/gridadaptinitializationindicator.hh rename to dumux/implicit/adaptive/gridadaptinitializationindicator.hh index cb63b396bb..ee992edde2 100644 --- a/dumux/implicit/common/gridadaptinitializationindicator.hh +++ b/dumux/implicit/adaptive/gridadaptinitializationindicator.hh @@ -25,12 +25,12 @@ /** * @file - * @brief Class defining an initialization indicator for grid adaption + * @brief Class defining an initialization indicator for grid adaptation */ namespace Dumux { /*!\ingroup ImplicitGridAdaptInitializationIndicator - * @brief Class defining an initialization indicator for grid adaption + * @brief Class defining an initialization indicator for grid adaptation * * Uses the defined grid adaptation indicator and further accounts for sources and boundaries. * Only for grid initialization! @@ -49,7 +49,7 @@ private: typedef typename GridView::Intersection Intersection; typedef typename GridView::template Codim<0>::Iterator ElementIterator; - typedef typename GET_PROP_TYPE(TypeTag, AdaptionIndicator) AdaptionIndicator; + typedef typename GET_PROP_TYPE(TypeTag, AdaptationIndicator) AdaptationIndicator; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; @@ -163,10 +163,10 @@ public: //First adapt for boundary conditions and sources to get a good initial solution if (nextMaxLevel_ == maxAllowedLevel_) - adaptionIndicator_.calculateIndicator(); + adaptationIndicator_.calculateIndicator(); // prepare an indicator for refinement - indicatorVector_.resize(problem_.model().numDofs()); + indicatorVector_.resize(problem_.gridView().size(0)); // set the default to coarsen indicatorVector_ = coarsenCell; @@ -241,7 +241,7 @@ public: if (indicatorVector_[idx] == refineCell) return true; else if (maxLevel_ == maxAllowedLevel_) - return adaptionIndicator_.refine(element); + return adaptationIndicator_.refine(element); else return false; } @@ -261,7 +261,7 @@ public: #endif if (indicatorVector_[idx] == coarsenCell && maxLevel_ < maxAllowedLevel_) return true; - else if (indicatorVector_[idx] == coarsenCell && !adaptionIndicator_.refine(element)) + else if (indicatorVector_[idx] == coarsenCell && !adaptationIndicator_.refine(element)) return true; else return false; @@ -272,7 +272,7 @@ public: return maxLevel_; } - /*! \brief Initializes the adaption indicator class */ + /*! \brief Initializes the adaptation indicator class */ void init() {}; @@ -281,7 +281,7 @@ public: return nextMaxLevel_ == maxAllowedLevel_; } - /*! \brief Constructs a GridAdaptionIndicator instance + /*! \brief Constructs a GridAdaptationIndicator instance * * This standard indicator is based on the saturation gradient. It checks the local gradient * compared to the maximum global gradient. The indicator is compared locally to a @@ -289,10 +289,10 @@ public: * or coarsening or should not be adapted. * * \param problem The problem object - * \param adaptionIndicator Indicator whether a be adapted + * \param adaptationIndicator Indicator whether a be adapted */ - ImplicitGridAdaptInitializationIndicator(Problem& problem, AdaptionIndicator& adaptionIndicator): - problem_(problem), adaptionIndicator_(adaptionIndicator), maxLevel_(0), nextMaxLevel_(0), eps_(1e-30) + ImplicitGridAdaptInitializationIndicator(Problem& problem, AdaptationIndicator& adaptationIndicator): + problem_(problem), adaptationIndicator_(adaptationIndicator), maxLevel_(0), nextMaxLevel_(0), eps_(1e-30) { minAllowedLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MinLevel); maxAllowedLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, GridAdapt, MaxLevel); @@ -310,7 +310,7 @@ public: private: Problem& problem_; - AdaptionIndicator& adaptionIndicator_; + AdaptationIndicator& adaptationIndicator_; Dune::DynamicVector<int> indicatorVector_; int maxLevel_; int nextMaxLevel_; @@ -325,7 +325,7 @@ private: /*!\ingroup IMPES - * @brief Class defining a start indicator for grid adaption + * @brief Class defining a start indicator for grid adaptation * *Default implementation * @@ -338,7 +338,7 @@ private: typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GridView::Traits::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, AdaptionIndicator) AdaptionIndicator; + typedef typename GET_PROP_TYPE(TypeTag, AdaptationIndicator) AdaptationIndicator; public: /*! \brief Calculates the indicator used for refinement/coarsening for each grid cell. @@ -374,18 +374,18 @@ public: return false; } - /*! \brief Initializes the adaption indicator class*/ + /*! \brief Initializes the adaptation indicator class*/ void init() {}; - /*! \brief Constructs a GridAdaptionIndicator for initialization of an adaptive grid + /*! \brief Constructs a GridAdaptationIndicator for initialization of an adaptive grid * * Default implementation * * \param problem The problem object - * \param adaptionIndicator Indicator whether a be adapted + * \param adaptationIndicator Indicator whether a be adapted */ - ImplicitGridAdaptInitializationIndicatorDefault(Problem& problem, AdaptionIndicator& adaptionIndicator) + ImplicitGridAdaptInitializationIndicatorDefault(Problem& problem, AdaptationIndicator& adaptationIndicator) {} }; diff --git a/dumux/implicit/common/gridadaptproperties.hh b/dumux/implicit/adaptive/gridadaptproperties.hh similarity index 91% rename from dumux/implicit/common/gridadaptproperties.hh rename to dumux/implicit/adaptive/gridadaptproperties.hh index fb0741d893..499316cf87 100644 --- a/dumux/implicit/common/gridadaptproperties.hh +++ b/dumux/implicit/adaptive/gridadaptproperties.hh @@ -33,19 +33,19 @@ namespace Dumux { namespace Properties { -//! Grid adaption type tag for all decoupled models. +//! Grid adaptation type tag for all decoupled models. NEW_TYPE_TAG(GridAdapt); //! Defines if the grid is h-adaptive NEW_PROP_TAG(AdaptiveGrid); //! Class defining the refinement/coarsening indicator -NEW_PROP_TAG(AdaptionIndicator); +NEW_PROP_TAG(AdaptationIndicator); //! Class defining the refinement/coarsening indicator for grid initialization -NEW_PROP_TAG(AdaptionInitializationIndicator); +NEW_PROP_TAG(AdaptationInitializationIndicator); -//! Switch the use of initial grid adaption on/off +//! Switch the use of initial grid adaptation on/off NEW_PROP_TAG(GridAdaptEnableInitializationIndicator); //! Mimimum allowed level @@ -66,8 +66,8 @@ NEW_PROP_TAG(GridAdaptRefineThreshold); //! Tolerance for coarsening NEW_PROP_TAG(GridAdaptCoarsenThreshold); -//! Time step interval for adaption -NEW_PROP_TAG(GridAdaptAdaptionInterval); +//! Time step interval for adaptation +NEW_PROP_TAG(GridAdaptAdaptationInterval); //! Switch for refinement at Dirichlet BC's -> not used by all indicators! NEW_PROP_TAG(GridAdaptRefineAtDirichletBC); diff --git a/dumux/implicit/common/gridadaptpropertydefaults.hh b/dumux/implicit/adaptive/gridadaptpropertydefaults.hh similarity index 87% rename from dumux/implicit/common/gridadaptpropertydefaults.hh rename to dumux/implicit/adaptive/gridadaptpropertydefaults.hh index 021814a64a..684c339de5 100644 --- a/dumux/implicit/common/gridadaptpropertydefaults.hh +++ b/dumux/implicit/adaptive/gridadaptpropertydefaults.hh @@ -47,8 +47,8 @@ SET_SCALAR_PROP(GridAdapt, GridAdaptRefineTolerance, 0.05); SET_SCALAR_PROP(GridAdapt, GridAdaptCoarsenTolerance, 0.001); SET_SCALAR_PROP(GridAdapt, GridAdaptRefineThreshold, 0.0); SET_SCALAR_PROP(GridAdapt, GridAdaptCoarsenThreshold, 0.0); -SET_INT_PROP(GridAdapt, GridAdaptAdaptionInterval, 1); -//Switch initial grid adaption off per default +SET_INT_PROP(GridAdapt, GridAdaptAdaptationInterval, 1); +//Switch initial grid adaptation off per default SET_BOOL_PROP(GridAdapt, GridAdaptEnableInitializationIndicator, false); // Switch of extra refinement strategy at boundaries/sources @@ -56,10 +56,10 @@ SET_BOOL_PROP(GridAdapt, GridAdaptRefineAtDirichletBC, false); SET_BOOL_PROP(GridAdapt, GridAdaptRefineAtFluxBC, false); SET_BOOL_PROP(GridAdapt, GridAdaptRefineAtSource, false); -//! Set the default indicator class models for adaption or coarsening -SET_TYPE_PROP(GridAdapt, AdaptionIndicator, ImplicitGridAdaptIndicatorDefault<TypeTag>); +//! Set the default indicator class models for adaptation or coarsening +SET_TYPE_PROP(GridAdapt, AdaptationIndicator, ImplicitGridAdaptIndicatorDefault<TypeTag>); //!Set default class for adaptation initialization indicator -SET_TYPE_PROP(GridAdapt, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicatorDefault<TypeTag>); +SET_TYPE_PROP(GridAdapt, AdaptationInitializationIndicator, ImplicitGridAdaptInitializationIndicatorDefault<TypeTag>); } // namespace Properties } // namespace Dumux diff --git a/dumux/implicit/common/implicitmodel.hh b/dumux/implicit/common/implicitmodel.hh index 5606585d37..098aa3b120 100644 --- a/dumux/implicit/common/implicitmodel.hh +++ b/dumux/implicit/common/implicitmodel.hh @@ -28,7 +28,7 @@ #include <dune/istl/bvector.hh> #include "implicitproperties.hh" -#include "gridadaptproperties.hh" +#include <dumux/implicit/adaptive/gridadaptproperties.hh> #include <dumux/common/valgrind.hh> #include <dumux/parallel/vertexhandles.hh> diff --git a/dumux/implicit/common/implicitproblem.hh b/dumux/implicit/common/implicitproblem.hh index 66eaeccddf..2d72c79c00 100644 --- a/dumux/implicit/common/implicitproblem.hh +++ b/dumux/implicit/common/implicitproblem.hh @@ -27,7 +27,7 @@ #include "implicitmodel.hh" #include <dumux/io/restart.hh> -#include <dumux/implicit/common/gridadapt.hh> +#include <dumux/implicit/adaptive/gridadapt.hh> namespace Dumux { @@ -920,7 +920,9 @@ private: // Tell the result writer that the grid changes if we are adaptive if (adaptiveGrid) + { resultWriter_->gridChanged(); + } } std::string simName_; @@ -945,6 +947,6 @@ private: }; } // namespace Dumux -#include <dumux/implicit/common/gridadaptpropertydefaults.hh> +#include <dumux/implicit/adaptive/gridadaptpropertydefaults.hh> #endif diff --git a/dumux/implicit/common/implicitproperties.hh b/dumux/implicit/common/implicitproperties.hh index 506f9ead69..4b8a7002c6 100644 --- a/dumux/implicit/common/implicitproperties.hh +++ b/dumux/implicit/common/implicitproperties.hh @@ -24,7 +24,7 @@ #include <dumux/common/basicproperties.hh> #include <dumux/linear/linearsolverproperties.hh> #include <dumux/nonlinear/newtonmethod.hh> -#include <dumux/implicit/common/gridadaptproperties.hh> +#include <dumux/implicit/adaptive/gridadaptproperties.hh> /*! * \ingroup Properties -- GitLab