From 5cbdfa7dda1d95551a6a84ec3fcd3a0f10eef4fa Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Tue, 6 Dec 2016 11:00:14 +0100 Subject: [PATCH 1/6] Perform the interface change of initialPhasePresence in a backward-compatible manner. Use SFINAE and enable_if to test if the old interface is still used in the problem file and to implement a corresponding wrapper function in `ImplicitProblem`. Once next is released, the wrappers and wrapper calls should be replaced by the new interface. --- dumux/discretization/box/subcontrolvolume.hh | 2 +- dumux/implicit/problem.hh | 78 +++++++++--- .../co2/implicit/volumevariables.hh | 34 ++--- .../compositional/primaryvariableswitch.hh | 4 +- .../co2/implicit/heterogeneousproblem.hh | 117 +++++++++--------- 5 files changed, 133 insertions(+), 102 deletions(-) diff --git a/dumux/discretization/box/subcontrolvolume.hh b/dumux/discretization/box/subcontrolvolume.hh index 00c9c29a1c..e65b840f08 100644 --- a/dumux/discretization/box/subcontrolvolume.hh +++ b/dumux/discretization/box/subcontrolvolume.hh @@ -85,7 +85,7 @@ public: return Geometry(Dune::GeometryType(Dune::GeometryType::cube, dim), corners_); } - //! The global index of this scv + //! The local index of this scv IndexType index() const { return scvIdx_; diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index 5bbeb00281..1c37e5ea0c 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -26,12 +26,52 @@ #include "properties.hh" #include "model.hh" +#include #include #include #include namespace Dumux { +// The following constructs are based on http://stackoverflow.com/a/16824239/6671978 +template +struct hasInitialPhasePresence {}; + +template +struct hasInitialPhasePresence { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().initialPhasePresence( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasInitialPhasePresenceAtPos {}; + +template +struct hasInitialPhasePresenceAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().initialPhasePresence( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + /*! * \ingroup ImplicitBaseProblems * \brief Base class for all fully implicit problems @@ -453,29 +493,31 @@ public: } /*! - * \brief Evaluate the initial phase state inside a control volume. + * \brief Returns the initial phase state for a sub control volume. * - * \param vertex The vertex - * \param vIdxGlobal The global index of the vertex - * \param globalPos The global position + * \param scv The sub control volume + * + * \todo replace by initialPhasePresence after release of next */ - int initialPhasePresence(const Vertex &vertex, - int &vIdxGlobal, - const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, int> + wrappedInitialPhasePresence(const SubControlVolume& scv) const { - // forward to generic interface - return asImp_().initialPhasePresenceAtPos(globalPos); + return asImp_().initialPhasePresence(scv); } - - /*! - * \brief Evaluate the initial value for a control volume. - * - * \param globalPos The global position - */ - void initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + template + DUNE_DEPRECATED_MSG("You are using the old interface 'initialPhasePresence(Vertex& ...)' in your problem file. Change to 'initialPhasePresence(SCV&)'.") + typename std::enable_if_t::value, int> + wrappedInitialPhasePresence(const SubControlVolume& scv) const { - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide a initialPhasePresenceAtPos() method."); + int dofIndex = scv.dofIndex(); + return asImp_().initialPhasePresence(*gridView_.template begin(), dofIndex, scv.dofPosition()); + } + template + typename std::enable_if_t::value, int> + wrappedInitialPhasePresence(const SubControlVolume& scv) const + { + return asImp_().initialPhasePresenceAtPos(scv.dofPosition()); } /*! diff --git a/dumux/porousmediumflow/co2/implicit/volumevariables.hh b/dumux/porousmediumflow/co2/implicit/volumevariables.hh index ff32876146..eead00c6af 100644 --- a/dumux/porousmediumflow/co2/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/co2/implicit/volumevariables.hh @@ -47,6 +47,7 @@ class CO2VolumeVariables: public TwoPTwoCVolumeVariables // typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; @@ -106,25 +107,14 @@ public: void update(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx, - const bool isOldSol) + const SubControlVolume& scv) { - // Update BoxVolVars but not 2p2cvolvars - // ToDo: Is BaseClassType the right name? - BaseClassType::update(priVars, - problem, - element, - fvGeometry, - scvIdx, - isOldSol); + BaseClassType::update(priVars, problem, element, scv); - int dofIdxGlobal = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim); + auto phasePresence = problem.model().priVarSwitch().phasePresence(scv.dofIndex()); - int phasePresence = problem.model().phasePresence(dofIdxGlobal, isOldSol); - - Scalar temp = Implementation::temperature_(priVars, problem, element, fvGeometry, scvIdx); - ParentType::fluidState_.setTemperature(temp); + Scalar t = Implementation::temperature(priVars, problem, element, scv); + ParentType::fluidState_.setTemperature(t); ///////////// // set the saturations @@ -148,7 +138,7 @@ public: // capillary pressure parameters const MaterialLawParams &materialParams = - problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx); + problem.spatialParams().materialLawParams(element, scv); Scalar pc = MaterialLaw::pc(materialParams, 1 - sn); @@ -323,9 +313,7 @@ public: } // porosity - ParentType::porosity_ = problem.spatialParams().porosity(element, - fvGeometry, - scvIdx); + ParentType::porosity_ = problem.spatialParams().porosity(element, scv); Valgrind::CheckDefined(ParentType::porosity_); // if(phasePresence == bothPhases) // { @@ -345,7 +333,7 @@ public: // } // energy related quantities not contained in the fluid state - asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol); + asImp_().updateEnergy_(priVars, problem, element, scv); } @@ -376,9 +364,7 @@ protected: void updateEnergy_(const PrimaryVariables &sol, const Problem &problem, const Element &element, - const FVElementGeometry &fvGeometry, - const int vIdx, - bool isOldSol) + const SubControlVolume& scv) { } diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh index e8b9a87d2f..a058fb22be 100644 --- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh +++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh @@ -25,6 +25,7 @@ #define DUMUX_PRIMARY_VARIABLE_SWITCH_HH #include +#include namespace Dumux { @@ -71,7 +72,8 @@ public: for (auto&& scv : scvs(fvGeometry)) { auto dofIdxGlobal = scv.dofIndex(); - phasePresence_[dofIdxGlobal] = problem.initialPhasePresence(scv); + // \todo replace by initialPhasePresence after release of next + phasePresence_[dofIdxGlobal] = problem.wrappedInitialPhasePresence(scv); } } diff --git a/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh b/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh index 155955e493..41233e9c51 100644 --- a/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh +++ b/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh @@ -24,6 +24,7 @@ #ifndef DUMUX_HETEROGENEOUS_PROBLEM_HH #define DUMUX_HETEROGENEOUS_PROBLEM_HH +#include #include #include #include @@ -44,7 +45,7 @@ namespace Properties { NEW_TYPE_TAG(HeterogeneousProblem, INHERITS_FROM(TwoPTwoC, HeterogeneousSpatialParams)); NEW_TYPE_TAG(HeterogeneousBoxProblem, INHERITS_FROM(BoxModel, HeterogeneousProblem)); -NEW_TYPE_TAG(HeterogeneousCCProblem, INHERITS_FROM(CCModel, HeterogeneousProblem)); +NEW_TYPE_TAG(HeterogeneousCCProblem, INHERITS_FROM(CCTpfaModel, HeterogeneousProblem)); // Set the grid type #if HAVE_DUNE_ALUGRID @@ -220,16 +221,16 @@ public: */ void postTimeStep() { - // Calculate storage terms - PrimaryVariables storageL, storageG; - this->model().globalPhaseStorage(storageL, lPhaseIdx); - this->model().globalPhaseStorage(storageG, gPhaseIdx); - - // Write mass balance information for rank 0 - if (this->gridView().comm().rank() == 0) { - std::cout<<"Storage: liquid=[" << storageL << "]" - << " gas=[" << storageG << "]\n"; - } +// // Calculate storage terms +// PrimaryVariables storageL, storageG; +// this->model().globalPhaseStorage(storageL, lPhaseIdx); +// this->model().globalPhaseStorage(storageG, gPhaseIdx); +// +// // Write mass balance information for rank 0 +// if (this->gridView().comm().rank() == 0) { +// std::cout<<"Storage: liquid=[" << storageL << "]" +// << " gas=[" << storageG << "]\n"; +// } } /*! @@ -237,53 +238,53 @@ public: * from the solution of the current time step to the VTK * writer. */ - void addOutputVtkFields() - { - typedef Dune::BlockVector > ScalarField; - - // get the number of degrees of freedom - unsigned numDofs = this->model().numDofs(); - unsigned numElements = this->gridView().size(0); - - //create required scalar fields - ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numElements); - ScalarField *cellPorosity = this->resultWriter().allocateManagedBuffer(numElements); - ScalarField *boxVolume = this->resultWriter().allocateManagedBuffer(numDofs); - (*boxVolume) = 0; - - //Fill the scalar fields with values - ScalarField *rank = this->resultWriter().allocateManagedBuffer(numElements); - - FVElementGeometry fvGeometry; - VolumeVariables volVars; - - for (const auto& element : elements(this->gridView())) - { - int eIdx = this->elementMapper().index(element); - (*rank)[eIdx] = this->gridView().comm().rank(); - fvGeometry.update(this->gridView(), element); - - - for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) - { - int dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim); - volVars.update(this->model().curSol()[dofIdxGlobal], - *this, - element, - fvGeometry, - scvIdx, - false); - (*boxVolume)[dofIdxGlobal] += fvGeometry.subContVol[scvIdx].volume; - } - (*Kxx)[eIdx] = this->spatialParams().intrinsicPermeability(element, fvGeometry, /*element data*/ 0); - (*cellPorosity)[eIdx] = this->spatialParams().porosity(element, fvGeometry, /*element data*/ 0); - } - - //pass the scalar fields to the vtkwriter - this->resultWriter().attachDofData(*Kxx, "Kxx", false); //element data - this->resultWriter().attachDofData(*cellPorosity, "cellwisePorosity", false); //element data - this->resultWriter().attachDofData(*boxVolume, "boxVolume", isBox); - } +// void addOutputVtkFields() +// { +// typedef Dune::BlockVector > ScalarField; +// +// // get the number of degrees of freedom +// unsigned numDofs = this->model().numDofs(); +// unsigned numElements = this->gridView().size(0); +// +// //create required scalar fields +// ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numElements); +// ScalarField *cellPorosity = this->resultWriter().allocateManagedBuffer(numElements); +// ScalarField *boxVolume = this->resultWriter().allocateManagedBuffer(numDofs); +// (*boxVolume) = 0; +// +// //Fill the scalar fields with values +// ScalarField *rank = this->resultWriter().allocateManagedBuffer(numElements); +// +// FVElementGeometry fvGeometry; +// VolumeVariables volVars; +// +// for (const auto& element : elements(this->gridView())) +// { +// int eIdx = this->elementMapper().index(element); +// (*rank)[eIdx] = this->gridView().comm().rank(); +// fvGeometry.update(this->gridView(), element); +// +// +// for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) +// { +// int dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim); +// volVars.update(this->model().curSol()[dofIdxGlobal], +// *this, +// element, +// fvGeometry, +// scvIdx, +// false); +// (*boxVolume)[dofIdxGlobal] += fvGeometry.subContVol[scvIdx].volume; +// } +// (*Kxx)[eIdx] = this->spatialParams().intrinsicPermeability(element, fvGeometry, /*element data*/ 0); +// (*cellPorosity)[eIdx] = this->spatialParams().porosity(element, fvGeometry, /*element data*/ 0); +// } +// +// //pass the scalar fields to the vtkwriter +// this->resultWriter().attachDofData(*Kxx, "Kxx", false); //element data +// this->resultWriter().attachDofData(*cellPorosity, "cellwisePorosity", false); //element data +// this->resultWriter().attachDofData(*boxVolume, "boxVolume", isBox); +// } /*! * \name Problem parameters -- GitLab From f7d0492a3365924060985b5f3f7a1d2a586095b8 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Wed, 7 Dec 2016 15:22:33 +0100 Subject: [PATCH 2/6] change name from `wrapped...` to `...Wrapper`, add `Element` to the parameter list --- dumux/implicit/problem.hh | 15 ++++++++------- .../compositional/primaryvariableswitch.hh | 2 +- .../2p2c/implicit/injectionproblem.hh | 2 +- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index 1c37e5ea0c..b42155a86f 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -495,27 +495,28 @@ public: /*! * \brief Returns the initial phase state for a sub control volume. * + * \param element The element * \param scv The sub control volume * * \todo replace by initialPhasePresence after release of next */ template - typename std::enable_if_t::value, int> - wrappedInitialPhasePresence(const SubControlVolume& scv) const + typename std::enable_if_t::value, int> + initialPhasePresenceWrapper(const Element& element, const SubControlVolume& scv) const { - return asImp_().initialPhasePresence(scv); + return asImp_().initialPhasePresence(element, scv); } template - DUNE_DEPRECATED_MSG("You are using the old interface 'initialPhasePresence(Vertex& ...)' in your problem file. Change to 'initialPhasePresence(SCV&)'.") + DUNE_DEPRECATED_MSG("You are using the old interface 'initialPhasePresence(Vertex& ...)' in your problem file. Change to 'initialPhasePresence(Element&, SCV&)'.") typename std::enable_if_t::value, int> - wrappedInitialPhasePresence(const SubControlVolume& scv) const + initialPhasePresenceWrapper(const Element& element, const SubControlVolume& scv) const { int dofIndex = scv.dofIndex(); - return asImp_().initialPhasePresence(*gridView_.template begin(), dofIndex, scv.dofPosition()); + return asImp_().initialPhasePresence(element.template subEntity(scv.index()), dofIndex, scv.dofPosition()); } template typename std::enable_if_t::value, int> - wrappedInitialPhasePresence(const SubControlVolume& scv) const + initialPhasePresenceWrapper(const Element& element, const SubControlVolume& scv) const { return asImp_().initialPhasePresenceAtPos(scv.dofPosition()); } diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh index a058fb22be..7bcd93c290 100644 --- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh +++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh @@ -73,7 +73,7 @@ public: { auto dofIdxGlobal = scv.dofIndex(); // \todo replace by initialPhasePresence after release of next - phasePresence_[dofIdxGlobal] = problem.wrappedInitialPhasePresence(scv); + phasePresence_[dofIdxGlobal] = problem.initialPhasePresenceWrapper(element, scv); } } diff --git a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh index 2335de5151..7300238467 100644 --- a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh +++ b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh @@ -330,7 +330,7 @@ public: * \param vIdxGlobal The global index of the vertex * \param globalPos The global position */ - int initialPhasePresence(const SubControlVolume& scv) const + int initialPhasePresence(const Element& element, const SubControlVolume& scv) const { return Indices::wPhaseOnly; } // \} -- GitLab From 09d44219bb9b0c0641872e41119ee478dc1233bb Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Thu, 8 Dec 2016 11:29:32 +0100 Subject: [PATCH 3/6] Perform the interface change of `initial` in a backward-compatible manner. Use SFINAE and enable_if to test if the old interface is still used in the problem file and to implement a corresponding wrapper function in `ImplicitProblem`. Once next is released, the wrappers and wrapper calls should be replaced by the new interface. If the old general interface `void initial(PrimaryVariables&, const Element&, const FVElementGeometry&, const int)` is used, this introduces a performance penalty. For every call, the corresponding `FVElementGeometry` is constructed. Conflicts: dumux/porousmediumflow/implicit/fluxvariables.hh --- dumux/implicit/model.hh | 2 +- dumux/implicit/problem.hh | 97 +++++++++++++------ .../co2/implicit/volumevariables.hh | 2 +- .../heterogeneousspatialparameters.hh | 11 +-- 4 files changed, 76 insertions(+), 36 deletions(-) diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index fb14955585..141d2a7679 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -812,7 +812,7 @@ protected: { // let the problem do the dirty work of nailing down // the initial solution. - auto initPriVars = problem_().initial(scv); + auto initPriVars = problem_().initialWrapper(element, scv); auto dofIdxGlobal = scv.dofIndex(); if (isBox) diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index b42155a86f..d1f9b6e38f 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -61,7 +61,45 @@ struct hasInitialPhasePresenceAtPos { private: template static constexpr auto check(T*) - -> typename std::is_same().initialPhasePresence( std::declval()... ) ), + -> typename std::is_same().initialPhasePresenceAtPos( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasInitial {}; + +template +struct hasInitial { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().initial( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasInitialAtPos {}; + +template +struct hasInitialAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().initialAtPos( std::declval()... ) ), Ret>::type; template static constexpr std::false_type check(...); @@ -457,39 +495,44 @@ public: const GlobalPosition &globalPos) const {} /*! - * \brief Evaluate the initial value for a control volume. + * \brief Returns the initial values for a sub control volume. * - * \param values The initial values for the primary variables - * \param element The finite element - * \param fvGeometry The finite-volume geometry - * \param scvIdx The local subcontrolvolume index + * \param element The element + * \param scv The sub control volume * - * For this method, the \a values parameter stores primary - * variables. + * \todo replace by initial after release of next */ - PrimaryVariables initial(const SubControlVolume &scv) const + template + typename std::enable_if_t::value, PrimaryVariables> + initialWrapper(const Element& element, const SubControlVolume& scv) const + { + return asImp_().initial(element, scv); + } + template + typename std::enable_if_t::value, PrimaryVariables> + initialWrapper(const Element& element, const SubControlVolume& scv) const { - // forward to generic interface return asImp_().initialAtPos(scv.dofPosition()); } - - /*! - * \brief Evaluate the initial value for a control volume. - * - * \param values The initial values for the primary variables - * \param globalPos The position of the center of the finite volume - * for which the initial values ought to be - * set (in global coordinates) - * - * For this method, the \a values parameter stores primary variables. - */ - PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void initial(...)' in your problem file. Change to 'PrimaryVariables initial(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + initialWrapper(const Element& element, const SubControlVolume& scv) const { - // Throw an exception (there is no reasonable default value - // for initial values) - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide " - "a initialAtPos() method."); + FVElementGeometry fvGeometry; + fvGeometry.update(gridView_, element); + PrimaryVariables values(0); + asImp_().initial(values, element, fvGeometry, scv.index()); + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void initialAtPos(...)' in your problem file. Change to 'PrimaryVariables initialAtPos(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + initialWrapper(const Element& element, const SubControlVolume& scv) const + { + PrimaryVariables values(0); + asImp_().initialAtPos(values, scv.dofPosition()); + return values; } /*! diff --git a/dumux/porousmediumflow/co2/implicit/volumevariables.hh b/dumux/porousmediumflow/co2/implicit/volumevariables.hh index eead00c6af..a217959c2f 100644 --- a/dumux/porousmediumflow/co2/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/co2/implicit/volumevariables.hh @@ -113,7 +113,7 @@ public: auto phasePresence = problem.model().priVarSwitch().phasePresence(scv.dofIndex()); - Scalar t = Implementation::temperature(priVars, problem, element, scv); + Scalar t = BaseClassType::temperature(priVars, problem, element, scv); ParentType::fluidState_.setTemperature(t); ///////////// diff --git a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh index 5f074bd09d..393856c395 100644 --- a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh +++ b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh @@ -82,6 +82,7 @@ class HeterogeneousSpatialParams : public ImplicitSpatialParams enum { dim=GridView::dimension }; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; + typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GridView::template Codim<0>::Entity Element; public: @@ -151,8 +152,7 @@ public: * \param scvIdx The local index of the sub-control volume */ const Scalar intrinsicPermeability(const Element &element, - const FVElementGeometry &fvGeometry, - int scvIdx) const + const SubControlVolume& scv) const { //Get the global index of the element int eIdx = gridView_.indexSet().index(element); @@ -172,9 +172,7 @@ public: * \param fvGeometry The finite volume geometry of the element * \param scvIdx The local index of the sub-control volume */ - Scalar porosity(const Element &element, - const FVElementGeometry &fvGeometry, - int scvIdx) const + Scalar porosity(const Element &element, const SubControlVolume& scv) const { //Get the global index of the element int eIdx = gridView_.indexSet().index(element); @@ -196,8 +194,7 @@ public: * \param scvIdx The local index of the sub-control volume */ const MaterialLawParams& materialLawParams(const Element &element, - const FVElementGeometry &fvGeometry, - int scvIdx) const + const SubControlVolume& scv) const { return materialParams_; } -- GitLab From 438adc1e825809dbff2bac616f252cc5eebdb2be Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Thu, 8 Dec 2016 14:18:52 +0100 Subject: [PATCH 4/6] Perform the interface change of `boundaryTypes` in a backward-compatible manner. Use SFINAE and enable_if to test if the old interface is still used in the problem file and to implement corresponding wrapper functions in `ImplicitProblem`. Once next is released, the wrappers and wrapper calls should be replaced by the new interface. If the old general interface `void boundaryTypes(BoundaryTypes&, const Intersection&)` is used, this introduces a performance penalty. For every call, the corresponding `Intersection` has to be found by iterating over all element intersections. Conflicts: dumux/implicit/cellcentered/localresidual.hh dumux/implicit/cellcentered/mpfa/localresidual.hh dumux/porousmediumflow/implicit/fluxvariables.hh --- .../cellcentered/globalvolumevariables.hh | 2 +- .../mpfa/globalvolumevariables.hh | 2 +- .../cellcentered/mpfa/helper.hh | 2 +- .../tpfa/elementvolumevariables.hh | 2 +- dumux/geomechanics/elastic/model.hh | 2 +- dumux/implicit/box/elementboundarytypes.hh | 2 +- dumux/implicit/box/intersectiontovertexbc.hh | 13 +- .../cellcentered/elementboundarytypes.hh | 2 +- .../cellcentered/elementvolumevariables.hh | 12 +- dumux/implicit/cellcentered/localresidual.hh | 4 +- .../cellcentered/mpfa/localresidual.hh | 4 +- .../cornerpoint/elementvolumevariables.hh | 12 +- dumux/implicit/problem.hh | 141 +++++++++++++----- .../multidomain/2cstokes2p2c/localoperator.hh | 2 +- dumux/porousmediumflow/2p2c/implicit/model.hh | 2 +- .../2pncmin/implicit/model.hh | 2 +- .../implicit/velocityoutput.hh | 34 +---- .../heterogeneousspatialparameters.hh | 5 +- 18 files changed, 143 insertions(+), 102 deletions(-) diff --git a/dumux/discretization/cellcentered/globalvolumevariables.hh b/dumux/discretization/cellcentered/globalvolumevariables.hh index a0b5bb91f2..85bd7cd30e 100644 --- a/dumux/discretization/cellcentered/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/globalvolumevariables.hh @@ -85,7 +85,7 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = problem.boundaryTypes(element, scvf); + const auto bcTypes = problem.boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { const auto insideScvIdx = scvf.insideScvIdx(); diff --git a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh index 0112162568..b85c52b164 100644 --- a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh @@ -86,7 +86,7 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = problem.boundaryTypes(element, scvf); + const auto bcTypes = problem.boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { const auto insideScvIdx = scvf.insideScvIdx(); diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index bc5a38b1ea..338fb0b19a 100644 --- a/dumux/discretization/cellcentered/mpfa/helper.hh +++ b/dumux/discretization/cellcentered/mpfa/helper.hh @@ -525,7 +525,7 @@ public: if (!scvf.boundary()) return MpfaFaceTypes::interior; - auto bcTypes = problem.boundaryTypes(element, scvf); + auto bcTypes = problem.boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyNeumann()) return MpfaFaceTypes::neumann; if (bcTypes.hasOnlyDirichlet()) diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index 168bb500fd..c8c6ae0d91 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -150,7 +150,7 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf); + const auto bcTypes = globalVolVars().problem_().boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf); diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh index b2c5c13b1b..18311bde1e 100644 --- a/dumux/geomechanics/elastic/model.hh +++ b/dumux/geomechanics/elastic/model.hh @@ -144,7 +144,7 @@ public: { if (scvFace.boundary()) { - BoundaryTypes bcTypes = this->problem_().boundaryTypes(element, scvFace); + BoundaryTypes bcTypes = this->problem_().boundaryTypesWrapper(element, scvFace); if (bcTypes.hasNeumann()) continue; } diff --git a/dumux/implicit/box/elementboundarytypes.hh b/dumux/implicit/box/elementboundarytypes.hh index 1a8e24d072..9cf36e28a2 100644 --- a/dumux/implicit/box/elementboundarytypes.hh +++ b/dumux/implicit/box/elementboundarytypes.hh @@ -96,7 +96,7 @@ public: if (problem.model().onBoundary(scv)) { - (*this)[scvIdxLocal] = problem.boundaryTypes(element, scv); + (*this)[scvIdxLocal] = problem.boundaryTypesWrapper(element, scv); hasDirichlet_ = hasDirichlet_ || (*this)[scvIdxLocal].hasDirichlet(); hasNeumann_ = hasNeumann_ || (*this)[scvIdxLocal].hasNeumann(); diff --git a/dumux/implicit/box/intersectiontovertexbc.hh b/dumux/implicit/box/intersectiontovertexbc.hh index 077a324193..da81a72076 100644 --- a/dumux/implicit/box/intersectiontovertexbc.hh +++ b/dumux/implicit/box/intersectiontovertexbc.hh @@ -38,6 +38,7 @@ class IntersectionToVertexBC typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; enum { @@ -57,22 +58,24 @@ public: vertexBC.resize(problem.vertexMapper().size()); for (int vIdx = 0; vIdx < vertexBC.size(); vIdx++) vertexBC[vIdx].setAllNeumann(); + const auto& globalFVGeometry = problem.model().globalFvGeometry(); + FVElementGeometry fvGeometry(globalFVGeometry); for (const auto& element : elements(problem.gridView())) { Dune::GeometryType geomType = element.geometry().type(); const ReferenceElement &refElement = ReferenceElements::general(geomType); + fvGeometry.bind(element); - for (const auto& intersection : intersections(problem.gridView(), element)) { - if (!intersection.boundary()) + for (const auto& scvf : scvfs(fvGeometry)) { + if (!scvf.boundary()) continue; - BoundaryTypes bcTypes; - problem.boundaryTypes(bcTypes, intersection); + auto bcTypes = problem.boundaryTypesWrapper(element, scvf); if (!bcTypes.hasDirichlet()) continue; - int fIdx = intersection.indexInInside(); + int fIdx = scvf.index(); int numFaceVerts = refElement.size(fIdx, 1, dim); for (int faceVertexIdx = 0; faceVertexIdx < numFaceVerts; ++faceVertexIdx) { diff --git a/dumux/implicit/cellcentered/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh index 606a9e503d..ef3a5f2a91 100644 --- a/dumux/implicit/cellcentered/elementboundarytypes.hh +++ b/dumux/implicit/cellcentered/elementboundarytypes.hh @@ -97,7 +97,7 @@ public: if (!scvFace.boundary()) continue; - (*this)[0] = problem.boundaryTypes(element, scvFace); + (*this)[0] = problem.boundaryTypesWrapper(element, scvFace); hasDirichlet_ = hasDirichlet_ || (*this)[0].hasDirichlet(); hasNeumann_ = hasNeumann_ || (*this)[0].hasNeumann(); diff --git a/dumux/implicit/cellcentered/elementvolumevariables.hh b/dumux/implicit/cellcentered/elementvolumevariables.hh index c09efacc7c..eb974e6dc4 100644 --- a/dumux/implicit/cellcentered/elementvolumevariables.hh +++ b/dumux/implicit/cellcentered/elementvolumevariables.hh @@ -87,19 +87,17 @@ public: isDirichlet_ = std::vector(element.subEntities(1), false); // add volume variables for the boundary faces - for (const auto& intersection : intersections(problem.gridView(), element)) + for (const auto& scvf : scvfs(fvGeometry)) { - if (!intersection.boundary()) + if (!scvf.boundary()) continue; - BoundaryTypes bcTypes; - problem.boundaryTypes(bcTypes, intersection); + auto bcTypes = problem.boundaryTypesWrapper(element, scvf); - int fIdx = intersection.indexInInside(); + int fIdx = scvf.index(); if (bcTypes.hasDirichlet()) { - PrimaryVariables dirichletValues; - problem.dirichlet(dirichletValues, intersection); + PrimaryVariables dirichletValues = problem.dirichlet(element, scvf); ghostVolVars_[fIdx].update(dirichletValues, problem, diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index 05d73ea4d9..3bbdade887 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -105,7 +105,7 @@ protected: { if (scvf.boundary()) { - auto bcTypes = this->problem().boundaryTypes(element, scvf); + auto bcTypes = this->problem().boundaryTypesWrapper(element, scvf); this->residual_[0] += evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache); } } @@ -158,7 +158,7 @@ protected: const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace &scvf) { - BoundaryTypes bcTypes = this->problem().boundaryTypes(element, scvf); + BoundaryTypes bcTypes = this->problem().boundaryTypesWrapper(element, scvf); if (bcTypes.hasDirichlet() && bcTypes.hasNeumann()) this->asImp_().evalDirichletSegmentMixed_(element, fvGeometry, elemVolVars, scvf, bcTypes); diff --git a/dumux/implicit/cellcentered/mpfa/localresidual.hh b/dumux/implicit/cellcentered/mpfa/localresidual.hh index 6fd5394927..a91ea46470 100644 --- a/dumux/implicit/cellcentered/mpfa/localresidual.hh +++ b/dumux/implicit/cellcentered/mpfa/localresidual.hh @@ -103,7 +103,7 @@ protected: return this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache); else { - auto bcTypes = this->problem().boundaryTypes(element, scvf); + auto bcTypes = this->problem().boundaryTypesWrapper(element, scvf); return this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache); } } @@ -121,7 +121,7 @@ protected: { if (scvf.boundary()) { - auto bcTypes = this->problem().boundaryTypes(element, scvf); + auto bcTypes = this->problem().boundaryTypesWrapper(element, scvf); this->residual_[0] += this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache); } } diff --git a/dumux/implicit/cornerpoint/elementvolumevariables.hh b/dumux/implicit/cornerpoint/elementvolumevariables.hh index 383200e7e7..5a23f3fce3 100644 --- a/dumux/implicit/cornerpoint/elementvolumevariables.hh +++ b/dumux/implicit/cornerpoint/elementvolumevariables.hh @@ -105,20 +105,18 @@ public: this->resize(numNeighbors + numFaces); // add volume variables for the boundary faces - for (const auto& intersection : intersections(problem.gridView(), element)) { - if (!intersection.boundary()) + for (const auto& scvf : scvfs(fvGeometry)) { + if (!scvf.boundary()) continue; - BoundaryTypes bcTypes; - problem.boundaryTypes(bcTypes, intersection); + auto bcTypes = problem.boundaryTypesWrapper(element, scvf); - int fIdx = intersection.indexInInside(); + int fIdx = scvf.index(); int indexInVariables = numNeighbors + fIdx; if (bcTypes.hasDirichlet()) { - PrimaryVariables dirichletValues; - problem.dirichlet(dirichletValues, intersection); + auto dirichletValues = problem.dirichlet(element, scvf); (*this)[indexInVariables].update(dirichletValues, problem, diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index d1f9b6e38f..bebf3c2996 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -110,6 +110,44 @@ public: static constexpr bool value = type::value; }; +template +struct hasBoundaryTypes {}; + +template +struct hasBoundaryTypes { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().boundaryTypes( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasBoundaryTypesAtPos {}; + +template +struct hasBoundaryTypesAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().boundaryTypesAtPos( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + /*! * \ingroup ImplicitBaseProblems * \brief Base class for all fully implicit problems @@ -228,55 +266,84 @@ public: } /*! - * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary segment. + * \brief Returns the boundary types for a sub control volume. + * + * \param element The element + * \param scv The sub control volume * - * \param values The boundary types for the conservation equations - * \param vertex The vertex for which the boundary type is set + * \todo replace by boundaryTypes after release of next */ - BoundaryTypes boundaryTypes(const Element &element, - const SubControlVolume &scv) const + template + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolume& scv) const + { + return asImp_().boundaryTypes(element, scv); + } + template + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolume& scv) const { - if (!isBox) - DUNE_THROW(Dune::InvalidStateException, - "boundaryTypes(..., scv) called for cell-centered method."); - - // forward it to the method which only takes the global coordinate return asImp_().boundaryTypesAtPos(scv.dofPosition()); } - - /*! - * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary segment. - * - * \param values The boundary types for the conservation equations - * \param intersection The intersection for which the boundary type is set - */ - BoundaryTypes boundaryTypes(const Element &element, - const SubControlVolumeFace &scvf) const + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void boundaryTypes(...)' in your problem file. Change to 'BoundaryTypes boundaryTypes(...)'.") + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolume& scv) const { - if (isBox) - DUNE_THROW(Dune::InvalidStateException, - "boundaryTypes(..., scvf) called for box method."); - - // forward it to the method which only takes the global coordinate - return asImp_().boundaryTypesAtPos(scvf.center()); + BoundaryTypes bcTypes; + asImp_().boundaryTypes(bcTypes, element.template subEntity(scv.index())); + return bcTypes; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void boundaryTypesAtPos(...)' in your problem file. Change to 'BoundaryTypes boundaryTypesAtPos(...)'.") + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolume& scv) const + { + BoundaryTypes bcTypes; + asImp_().boundaryTypesAtPos(bcTypes, scv.dofPosition()); + return bcTypes; } /*! - * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary segment. + * \brief Returns the boundary types for a sub control volume face. + * + * \param element The element + * \param scvf The sub control volume face * - * \param values The boundary types for the conservation equations - * \param globalPos The position of the finite volume in global coordinates + * \todo replace by boundaryTypes after release of next */ - BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolumeFace& scvf) const { - // Throw an exception (there is no reasonable default value - // for Dirichlet conditions) - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide " - "a boundaryTypes() method."); + return asImp_().boundaryTypes(element, scvf); + } + template + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + return asImp_().boundaryTypesAtPos(scvf.ipGlobal()); + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void boundaryTypes(...)' in your problem file. Change to 'BoundaryTypes boundaryTypes(...)'.") + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + BoundaryTypes bcTypes; + // The following only works for conforming grids: + for (const auto& intersection : intersections(gridView_, element)) + if (intersection.indexInInside() == scvf.index()) + asImp_().boundaryTypes(bcTypes, intersection); + return bcTypes; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void boundaryTypesAtPos(...)' in your problem file. Change to 'BoundaryTypes boundaryTypesAtPos(...)'.") + typename std::enable_if_t::value, BoundaryTypes> + boundaryTypesWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + BoundaryTypes bcTypes; + asImp_().boundaryTypesAtPos(bcTypes, scvf.ipGlobal()); + return bcTypes; } /*! diff --git a/dumux/multidomain/2cstokes2p2c/localoperator.hh b/dumux/multidomain/2cstokes2p2c/localoperator.hh index 32c54e2e7a..bde8cc0e22 100644 --- a/dumux/multidomain/2cstokes2p2c/localoperator.hh +++ b/dumux/multidomain/2cstokes2p2c/localoperator.hh @@ -299,7 +299,7 @@ public: const VertexPointer1 vPtr1 = (*sdElement1).template subEntity(vertInElem1); const VertexPointer2 vPtr2 = (*sdElement2).template subEntity(vertInElem2); - globalProblem_.sdProblem1().boundaryTypes(cParams.boundaryTypes1, vPtr1); + cParams.boundaryTypes1globalProblem_.sdProblem1().boundaryTypes(cParams.boundaryTypes1, vPtr1); globalProblem_.sdProblem2().boundaryTypes(cParams.boundaryTypes2, vPtr2); BoundaryVariables1 boundaryVars1; diff --git a/dumux/porousmediumflow/2p2c/implicit/model.hh b/dumux/porousmediumflow/2p2c/implicit/model.hh index 02c151390e..4160c36a32 100644 --- a/dumux/porousmediumflow/2p2c/implicit/model.hh +++ b/dumux/porousmediumflow/2p2c/implicit/model.hh @@ -188,7 +188,7 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = this->problem_().boundaryTypes(element, scvf); + const auto bcTypes = this->problem_().boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { const auto insideScvIdx = scvf.insideScvIdx(); diff --git a/dumux/porousmediumflow/2pncmin/implicit/model.hh b/dumux/porousmediumflow/2pncmin/implicit/model.hh index 650a5f10fc..de805f723e 100644 --- a/dumux/porousmediumflow/2pncmin/implicit/model.hh +++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh @@ -195,7 +195,7 @@ public: continue; // check if boundary is a pure dirichlet boundary - const auto bcTypes = this->problem_().boundaryTypes(element, scvf); + const auto bcTypes = this->problem_().boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { const auto insideScvIdx = scvf.insideScvIdx(); diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh index 8f386f1042..71c8ed9728 100644 --- a/dumux/porousmediumflow/implicit/velocityoutput.hh +++ b/dumux/porousmediumflow/implicit/velocityoutput.hh @@ -108,31 +108,6 @@ public: return velocityOutput_; } - // The following SFINAE enable_if usage allows compilation, even if only a - // - // boundaryTypes(BoundaryTypes&, const Vertex&) - // - // is provided in the problem file. In that case, the compiler cannot detect - // (without additional measures like "using...") the signature - // - // boundaryTypes(BoundaryTypes&, const Intersection&) - // - // in the problem base class. Therefore, calls to this method trigger a - // compiler error. However, that call is needed for calculating velocities - // if the cell-centered discretization is used. By proceeding as in the - // following lines, that call will only be compiled if cell-centered - // actually is used. - template - void problemBoundaryTypes(BoundaryTypes& bcTypes, - const typename std::enable_if::type& intersection) const - { - problem_.boundaryTypes(bcTypes, intersection); - } - template - void problemBoundaryTypes(BoundaryTypes& bcTypes, - const typename std::enable_if::type& intersection) const - {} - template void calculateVelocity(VelocityVector& velocity, const ElementVolumeVariables& elemVolVars, @@ -255,16 +230,15 @@ public: // Neumann conditions. if (element.hasBoundaryIntersections()) { - for (const auto& intersection : intersections(problem_.gridView(), element)) + for (const auto& scvf : scvfs(fvGeometry)) { - if (intersection.boundary()) + if (scvf.boundary()) { - BoundaryTypes bcTypes; - problemBoundaryTypes(bcTypes, intersection); + auto bcTypes = boundaryTypesWrapper(element, scvf); if (bcTypes.hasNeumann()) { - int fIdx = intersection.indexInInside(); + int fIdx = scvf.index(); // cubes if (dim == 1 || geomType.isCube()){ diff --git a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh index 393856c395..4a54910984 100644 --- a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh +++ b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh @@ -73,6 +73,7 @@ template class HeterogeneousSpatialParams : public ImplicitSpatialParams { typedef ImplicitSpatialParams ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; @@ -94,8 +95,8 @@ public: * * \param gridView The grid view */ - HeterogeneousSpatialParams(const GridView &gridView) - : ParentType(gridView), gridView_(gridView) + HeterogeneousSpatialParams(const Problem& problem, const GridView &gridView) + : ParentType(problem, gridView), gridView_(gridView) { /* * Layer Index Setup: -- GitLab From 4e841faf952ff345e805628aca3c9f625ae7ea81 Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Thu, 8 Dec 2016 16:42:24 +0100 Subject: [PATCH 5/6] Perform the interface change of `dirichlet` and `neumann` in a backward-compatible manner. Analogous comments as the ones for 714c973f and c8decaa5 apply. --- .../cellcentered/globalvolumevariables.hh | 2 +- .../mpfa/elementvolumevariables.hh | 4 +- .../mpfa/globalvolumevariables.hh | 2 +- .../mpfa/omethod/interactionvolume.hh | 2 +- .../tpfa/elementvolumevariables.hh | 2 +- dumux/implicit/box/localresidual.hh | 4 +- .../cellcentered/elementvolumevariables.hh | 2 +- dumux/implicit/cellcentered/localresidual.hh | 4 +- .../cornerpoint/elementvolumevariables.hh | 2 +- dumux/implicit/localresidual.hh | 2 +- dumux/implicit/problem.hh | 433 ++++++++++++++---- dumux/porousmediumflow/2p2c/implicit/model.hh | 2 +- .../2pncmin/implicit/model.hh | 2 +- .../mpnc/implicit/localresidual.hh | 2 +- .../heterogeneousspatialparameters.hh | 7 +- 15 files changed, 360 insertions(+), 112 deletions(-) diff --git a/dumux/discretization/cellcentered/globalvolumevariables.hh b/dumux/discretization/cellcentered/globalvolumevariables.hh index 85bd7cd30e..e974b12228 100644 --- a/dumux/discretization/cellcentered/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/globalvolumevariables.hh @@ -90,7 +90,7 @@ public: { const auto insideScvIdx = scvf.insideScvIdx(); const auto& insideScv = fvGeometry.scv(insideScvIdx); - const auto dirichletPriVars = problem.dirichlet(element, scvf); + const auto dirichletPriVars = problem.dirichletWrapper(element, scvf); volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv); } diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 03e8db6eb0..03ca9919c1 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -177,7 +177,7 @@ public: { // boundary volume variables VolumeVariables dirichletVolVars; - const auto dirichletPriVars = problem.dirichlet(element, scvf); + const auto dirichletPriVars = problem.dirichletWrapper(element, scvf); dirichletVolVars.update(dirichletPriVars, problem, element, scvI); volumeVariables_.emplace_back(std::move(dirichletVolVars)); @@ -215,7 +215,7 @@ public: // boundary volume variables VolumeVariables dirichletVolVars; auto&& ivScv = fvGeometry.scv(insideScvIdx); - const auto dirichletPriVars = problem.dirichlet(ivElement, ivScvf); + const auto dirichletPriVars = problem.dirichletWrapper(ivElement, ivScvf); dirichletVolVars.update(dirichletPriVars, problem, ivElement, ivScv); volumeVariables_.emplace_back(std::move(dirichletVolVars)); diff --git a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh index b85c52b164..0aa528a2bf 100644 --- a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh @@ -91,7 +91,7 @@ public: { const auto insideScvIdx = scvf.insideScvIdx(); const auto& insideScv = fvGeometry.scv(insideScvIdx); - const auto dirichletPriVars = problem.dirichlet(element, scvf); + const auto dirichletPriVars = problem.dirichletWrapper(element, scvf); volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv); } diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index fc476e6f90..d5cdedc858 100644 --- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh +++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh @@ -218,7 +218,7 @@ public: { const auto& element = localElement_(localScvf.insideLocalScvIndex()); const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex()); - auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx]; + auto neumannFlux = problem_().neumannWrapper(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx]; neumannFlux *= globalScvf.area(); // recover -k*gradh diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index c8c6ae0d91..2d12a92710 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -153,7 +153,7 @@ public: const auto bcTypes = globalVolVars().problem_().boundaryTypesWrapper(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf); + const auto dirichletPriVars = globalVolVars().problem_().dirichletWrapper(element, scvf); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); diff --git a/dumux/implicit/box/localresidual.hh b/dumux/implicit/box/localresidual.hh index 3dcb322ecb..bac97a199b 100644 --- a/dumux/implicit/box/localresidual.hh +++ b/dumux/implicit/box/localresidual.hh @@ -135,7 +135,7 @@ protected: const SubControlVolume &scv, const BoundaryTypes& bcTypes) { - PrimaryVariables dirichletValues = this->problem().dirichlet(element, scv); + PrimaryVariables dirichletValues = this->problem().dirichletWrapper(element, scv); // set the dirichlet conditions for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) @@ -188,7 +188,7 @@ protected: // temporary vector to store the neumann boundary fluxes PrimaryVariables flux(0); - auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf); + auto neumannFluxes = this->problem().neumannWrapper(element, fvGeometry, elemVolVars, scvf); // multiply neumann fluxes with the area and the extrusion factor neumannFluxes *= scvf.area()*elemVolVars[insideScv].extrusionFactor(); diff --git a/dumux/implicit/cellcentered/elementvolumevariables.hh b/dumux/implicit/cellcentered/elementvolumevariables.hh index eb974e6dc4..719dc8a9e5 100644 --- a/dumux/implicit/cellcentered/elementvolumevariables.hh +++ b/dumux/implicit/cellcentered/elementvolumevariables.hh @@ -97,7 +97,7 @@ public: int fIdx = scvf.index(); if (bcTypes.hasDirichlet()) { - PrimaryVariables dirichletValues = problem.dirichlet(element, scvf); + PrimaryVariables dirichletValues = problem.dirichletWrapper(element, scvf); ghostVolVars_[fIdx].update(dirichletValues, problem, diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index 3bbdade887..3b2b386e95 100644 --- a/dumux/implicit/cellcentered/localresidual.hh +++ b/dumux/implicit/cellcentered/localresidual.hh @@ -176,7 +176,7 @@ protected: // temporary vector to store the neumann boundary fluxes PrimaryVariables flux(0); - auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf); + auto neumannFluxes = this->problem().neumannWrapper(element, fvGeometry, elemVolVars, scvf); // multiply neumann fluxes with the area and the extrusion factor auto&& scv = fvGeometry.scv(scvf.insideScvIdx()); @@ -225,7 +225,7 @@ protected: const BoundaryTypes &bcTypes) { // temporary vector to store the Dirichlet values - PrimaryVariables dirichletValues = this->problem().dirichlet(element, scvf); + PrimaryVariables dirichletValues = this->problem().dirichletWrapper(element, scvf); // get the primary variables const auto& priVars = elemVolVars[scvf.insideScvIdx()].priVars(); diff --git a/dumux/implicit/cornerpoint/elementvolumevariables.hh b/dumux/implicit/cornerpoint/elementvolumevariables.hh index 5a23f3fce3..c3508d1e38 100644 --- a/dumux/implicit/cornerpoint/elementvolumevariables.hh +++ b/dumux/implicit/cornerpoint/elementvolumevariables.hh @@ -116,7 +116,7 @@ public: if (bcTypes.hasDirichlet()) { - auto dirichletValues = problem.dirichlet(element, scvf); + auto dirichletValues = problem.dirichletWrapper(element, scvf); (*this)[indexInVariables].update(dirichletValues, problem, diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index ee68de2417..b62f598634 100644 --- a/dumux/implicit/localresidual.hh +++ b/dumux/implicit/localresidual.hh @@ -230,7 +230,7 @@ public: PrimaryVariables source(0); // add contributions from volume flux sources - source += this->problem().source(element, fvGeometry, elemVolVars, scv); + source += this->problem().sourceWrapper(element, fvGeometry, elemVolVars, scv); // add contribution from possible point sources source += this->problem().scvPointSources(element, fvGeometry, elemVolVars, scv); diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh index bebf3c2996..89cab6bb86 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -148,6 +148,158 @@ public: static constexpr bool value = type::value; }; +template +struct hasDirichlet {}; + +template +struct hasDirichlet { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().dirichlet( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasDirichletAtPos {}; + +template +struct hasDirichletAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().dirichletAtPos( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasNeumann {}; + +template +struct hasNeumann { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().neumann( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasSolDependentNeumann {}; + +template +struct hasSolDependentNeumann { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().solDependentNeumann( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasNeumannAtPos {}; + +template +struct hasNeumannAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().neumannAtPos( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasSource {}; + +template +struct hasSource { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().source( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasSolDependentSource {}; + +template +struct hasSolDependentSource { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().solDependentSource( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + +template +struct hasSourceAtPos {}; + +template +struct hasSourceAtPos { +private: + template + static constexpr auto check(T*) + -> typename std::is_same().sourceAtPos( std::declval()... ) ), + Ret>::type; + + template static constexpr std::false_type check(...); + + typedef decltype(check(0)) type; + +public: + static constexpr bool value = type::value; +}; + /*! * \ingroup ImplicitBaseProblems * \brief Base class for all fully implicit problems @@ -333,7 +485,10 @@ public: // The following only works for conforming grids: for (const auto& intersection : intersections(gridView_, element)) if (intersection.indexInInside() == scvf.index()) + { asImp_().boundaryTypes(bcTypes, intersection); + break; + } return bcTypes; } template @@ -346,55 +501,89 @@ public: return bcTypes; } + /*! - * \brief Evaluate the boundary conditions for a dirichlet - * control volume. + * \brief Returns the Dirichlet values for a sub control volume. * - * \param values The dirichlet values for the primary variables - * \param scvFace the sub control volume face + * \param element The element + * \param scv The sub control volume * - * The method returns the boundary types information. + * \todo replace by dirichlet after release of next */ - PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + template + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolume& scv) const { - // forward it to the method which only takes the global coordinate - if (isBox) - { - DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method."); - } - else - return asImp_().dirichletAtPos(scvf.center()); + return asImp_().dirichlet(element, scv); } - - PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const + template + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolume& scv) const { - // forward it to the method which only takes the global coordinate - if (!isBox) - { - DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for cell-centered method."); - } - else - return asImp_().dirichletAtPos(scv.dofPosition()); + return asImp_().dirichletAtPos(scv.dofPosition()); + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void dirichlet(...)' in your problem file. Change to 'PrimaryVariables dirichlet(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolume& scv) const + { + PrimaryVariables values; + asImp_().dirichlet(values, element.template subEntity(scv.index())); + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void dirichletAtPos(...)' in your problem file. Change to 'PrimaryVariables dirichletAtPos(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolume& scv) const + { + PrimaryVariables values; + asImp_().dirichletAtPos(values, scv.dofPosition()); + return values; } /*! - * \brief Evaluate the boundary conditions for a dirichlet - * control volume. + * \brief Returns the Dirichlet values for a sub control volume face. * - * \param globalPos The position of the center of the finite volume - * for which the dirichlet condition ought to be - * set in global coordinates + * \param element The element + * \param scvf The sub control volume face * - * For this method, the \a values parameter stores primary variables. + * \todo replace by dirichlet after release of next */ - PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + return asImp_().dirichlet(element, scvf); + } + template + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolumeFace& scvf) const { - // Throw an exception (there is no reasonable default value - // for Dirichlet conditions) - DUNE_THROW(Dune::InvalidStateException, - "The problem specifies that some boundary " - "segments are dirichlet, but does not provide " - "a dirichlet() method."); + return asImp_().dirichletAtPos(scvf.ipGlobal()); + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void dirichlet(...)' in your problem file. Change to 'PrimaryVariables dirichlet(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + PrimaryVariables values; + // The following only works for conforming grids: + for (const auto& intersection : intersections(gridView_, element)) + if (intersection.indexInInside() == scvf.index()) + { + asImp_().dirichlet(values, intersection); + break; + } + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void dirichletAtPos(...)' in your problem file. Change to 'PrimaryVariables dirichletAtPos(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + dirichletWrapper(const Element& element, const SubControlVolumeFace& scvf) const + { + PrimaryVariables values; + asImp_().dirichletAtPos(values, scvf.ipGlobal()); + return values; } /*! @@ -405,48 +594,82 @@ public: * potentially solution dependent and requires some quantities that * are specific to the fully-implicit method. * - * \param values The neumann values for the conservation equations in units of + * \returns The neumann values for the conservation equations in units of * \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$ * \param element The finite element * \param fvGeometry The finite-volume geometry - * \param intersection The intersection between element and boundary - * \param scvIdx The local subcontrolvolume index - * \param boundaryFaceIdx The index of the boundary face * \param elemVolVars All volume variables for the element + * \param scvf The sub control volume face * * For this method, the \a values parameter stores the flux * in normal direction of each phase. Negative values mean influx. * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$. */ - PrimaryVariables neumann(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolvars, - const SubControlVolumeFace& scvFace) const + template + typename std::enable_if_t::value, PrimaryVariables> + neumannWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const { - // forward it to the interface with only the global position - return asImp_().neumannAtPos(scvFace.center()); + return asImp_().neumann(element, fvGeometry, elemVolVars, scvf); } - - /*! - * \brief Evaluate the boundary conditions for a neumann - * boundary segment. - * - * \param values The neumann values for the conservation equations in units of - * \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$ - * \param globalPos The position of the boundary face's integration point in global coordinates - * - * For this method, the \a values parameter stores the flux - * in normal direction of each phase. Negative values mean influx. - * E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$. - */ - PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, PrimaryVariables> + neumannWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const { - // Throw an exception (there is no reasonable default value - // for Neumann conditions) - DUNE_THROW(Dune::InvalidStateException, - "The problem specifies that some boundary " - "segments are neumann, but does not provide " - "a neumannAtPos() method."); + return asImp_().neumannAtPos(scvf.ipGlobal()); + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void solDependentNeumann(...)' in your problem file. Change to 'PrimaryVariables neumann(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + neumannWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const + { + PrimaryVariables values(0); + int scvIdx = 0; + for (const auto& intersection : intersections(gridView_, element)) + if (intersection.indexInInside() == scvf.index()) + { + asImp_().solDependentNeumann(values, element, fvGeometry, intersection, scvIdx, scvf.index(), elemVolVars); + break; + } + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void neumann(...)' in your problem file. Change to 'PrimaryVariables neumann(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + neumannWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const + { + PrimaryVariables values(0); + int scvIdx = 0; + for (const auto& intersection : intersections(gridView_, element)) + if (intersection.indexInInside() == scvf.index()) + { + asImp_().neumann(values, element, fvGeometry, intersection, scvIdx, scvf.index()); + break; + } + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void neumannAtPos(...)' in your problem file. Change to 'PrimaryVariables neumannAtPos(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + neumannWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const + { + PrimaryVariables values(0); + asImp_().neumannAtPos(values, scvf.ipGlobal()); + return values; } /*! @@ -457,7 +680,7 @@ public: * potentially solution dependent and requires some quantities that * are specific to the fully-implicit method. * - * \param values The source and sink values for the conservation equations in units of + * \returns The source and sink values for the conservation equations in units of * \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$ * \param element The finite element * \param fvGeometry The finite-volume geometry @@ -469,35 +692,59 @@ public: * that the conserved quantity is created, negative ones mean that it vanishes. * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. */ - PrimaryVariables source(const Element &element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolume &scv) const + template + typename std::enable_if_t::value, PrimaryVariables> + sourceWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume& scv) const + { + return asImp_().source(element, fvGeometry, elemVolVars, scv); + } + template + typename std::enable_if_t::value, PrimaryVariables> + sourceWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume& scv) const { - // forward to solution independent, fully-implicit specific interface return asImp_().sourceAtPos(scv.dofPosition()); } - - /*! - * \brief Evaluate the source term for all phases within a given - * sub-control-volume. - * - * \param values The source and sink values for the conservation equations in units of - * \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$ - * \param globalPos The position of the center of the finite volume - * for which the source term ought to be - * specified in global coordinates - * - * For this method, the \a values parameter stores the conserved quantity rate - * generated or annihilate per volume unit. Positive values mean - * that the conserved quantity is created, negative ones mean that it vanishes. - * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. - */ - PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void solDependentSource(...)' in your problem file. Change to 'PrimaryVariables source(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + sourceWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume& scv) const + { + PrimaryVariables values(0); + asImp_().solDependentSource(values, element, fvGeometry, scv.index(), elemVolVars); + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void source(...)' in your problem file. Change to 'PrimaryVariables source(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + sourceWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume& scv) const { - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide " - "a sourceAtPos() method."); + PrimaryVariables values(0); + asImp_().source(values, element, fvGeometry, scv.index()); + return values; + } + template + DUNE_DEPRECATED_MSG("You are using the old interface 'void sourceAtPos(...)' in your problem file. Change to 'PrimaryVariables sourceAtPos(...)'.") + typename std::enable_if_t::value, PrimaryVariables> + sourceWrapper(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume& scv) const + { + PrimaryVariables values(0); + asImp_().sourceAtPos(values, scv.dofPosition()); + return values; } /*! @@ -586,8 +833,8 @@ public: typename std::enable_if_t::value, PrimaryVariables> initialWrapper(const Element& element, const SubControlVolume& scv) const { - FVElementGeometry fvGeometry; - fvGeometry.update(gridView_, element); + FVElementGeometry fvGeometry(model_.globalFvGeometry()); + fvGeometry.bind(element); PrimaryVariables values(0); asImp_().initial(values, element, fvGeometry, scv.index()); return values; diff --git a/dumux/porousmediumflow/2p2c/implicit/model.hh b/dumux/porousmediumflow/2p2c/implicit/model.hh index 4160c36a32..b3162ed8ec 100644 --- a/dumux/porousmediumflow/2p2c/implicit/model.hh +++ b/dumux/porousmediumflow/2p2c/implicit/model.hh @@ -193,7 +193,7 @@ public: { const auto insideScvIdx = scvf.insideScvIdx(); const auto& insideScv = fvGeometry.scv(insideScvIdx); - const auto dirichletPriVars = this->problem_().dirichlet(element, scvf); + const auto dirichletPriVars = this->problem_().dirichletWrapper(element, scvf); this->nonConstCurGlobalVolVars().volVars(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), element, insideScv); } diff --git a/dumux/porousmediumflow/2pncmin/implicit/model.hh b/dumux/porousmediumflow/2pncmin/implicit/model.hh index de805f723e..7fa8131f0e 100644 --- a/dumux/porousmediumflow/2pncmin/implicit/model.hh +++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh @@ -200,7 +200,7 @@ public: { const auto insideScvIdx = scvf.insideScvIdx(); const auto& insideScv = fvGeometry.scv(insideScvIdx); - const auto dirichletPriVars = this->problem_().dirichlet(element, scvf); + const auto dirichletPriVars = this->problem_().dirichletWrapper(element, scvf); this->nonConstCurGlobalVolVars().volVars(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), element, insideScv); } diff --git a/dumux/porousmediumflow/mpnc/implicit/localresidual.hh b/dumux/porousmediumflow/mpnc/implicit/localresidual.hh index 3e7fd6e8ca..4ff2eb94b0 100644 --- a/dumux/porousmediumflow/mpnc/implicit/localresidual.hh +++ b/dumux/porousmediumflow/mpnc/implicit/localresidual.hh @@ -280,7 +280,7 @@ public: // temporary vector to store the Dirichlet boundary fluxes PrimaryVariables values; Valgrind::SetUndefined(values); - this->problem_().dirichlet(values, *isIt); + this->problem_().dirichletWrapper(values, *isIt); Valgrind::CheckDefined(values); // set Dirichlet conditions in a strong sense diff --git a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh index 4a54910984..8448788035 100644 --- a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh +++ b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh @@ -78,6 +78,7 @@ class HeterogeneousSpatialParams : public ImplicitSpatialParams typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; enum { dim=GridView::dimension }; @@ -152,11 +153,11 @@ public: * \param fvGeometry The finite volume geometry of the element * \param scvIdx The local index of the sub-control volume */ - const Scalar intrinsicPermeability(const Element &element, - const SubControlVolume& scv) const + Scalar intrinsicPermeability(const SubControlVolume& scv, + const VolumeVariables& volVars) const { //Get the global index of the element - int eIdx = gridView_.indexSet().index(element); + int eIdx = scv.elementIndex(); if (paramIdx_[eIdx] == barrierTop_) return barrierTopK_; -- GitLab From 92a759cc59f15894f961c1f2bd55a088c9e52e9b Mon Sep 17 00:00:00 2001 From: Bernd Flemisch Date: Fri, 9 Dec 2016 08:59:31 +0100 Subject: [PATCH 6/6] fix typo introduced in 714c973f --- dumux/multidomain/2cstokes2p2c/localoperator.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/multidomain/2cstokes2p2c/localoperator.hh b/dumux/multidomain/2cstokes2p2c/localoperator.hh index bde8cc0e22..32c54e2e7a 100644 --- a/dumux/multidomain/2cstokes2p2c/localoperator.hh +++ b/dumux/multidomain/2cstokes2p2c/localoperator.hh @@ -299,7 +299,7 @@ public: const VertexPointer1 vPtr1 = (*sdElement1).template subEntity(vertInElem1); const VertexPointer2 vPtr2 = (*sdElement2).template subEntity(vertInElem2); - cParams.boundaryTypes1globalProblem_.sdProblem1().boundaryTypes(cParams.boundaryTypes1, vPtr1); + globalProblem_.sdProblem1().boundaryTypes(cParams.boundaryTypes1, vPtr1); globalProblem_.sdProblem2().boundaryTypes(cParams.boundaryTypes2, vPtr2); BoundaryVariables1 boundaryVars1; -- GitLab