diff --git a/dumux/discretization/box/subcontrolvolume.hh b/dumux/discretization/box/subcontrolvolume.hh index 00c9c29a1c2dd8242c79c52fe4538036b4e767b2..e65b840f081b8b044e08239fc5b356d0a84f3f3e 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/discretization/cellcentered/globalvolumevariables.hh b/dumux/discretization/cellcentered/globalvolumevariables.hh index a0b5bb91f279323d9b5fa241a3d8c7381fa2b6a7..e974b12228521b4b6d414c18a1feb63cc6015a46 100644 --- a/dumux/discretization/cellcentered/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/globalvolumevariables.hh @@ -85,12 +85,12 @@ 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(); 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 03e8db6eb02e8aa80267f096644c53874adf8b4b..03ca9919c12e72cfac5bff382b2573c9bb8ff6e3 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 01121625683b66878a5d219d807cc0c0dc5addf4..0aa528a2bffd380323614d02c504044511a36dba 100644 --- a/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/globalvolumevariables.hh @@ -86,12 +86,12 @@ 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(); 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/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh index bc5a38b1ea486a452c70c7219ce86df30add9ef4..338fb0b19ac929d838dda54b07f63970865cd350 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/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh index fc476e6f90bf29e9b9d41271411a38974d22d1e2..d5cdedc858438e68c43ef077f62731dd8b3cd5db 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 168bb500fd80db32ec3b581d3ef7fd487ea99bae..2d12a9271050354ccf853cd22c847a169f780677 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -150,10 +150,10 @@ 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); + const auto dirichletPriVars = globalVolVars().problem_().dirichletWrapper(element, scvf); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh index b2c5c13b1bbe0cccaabc0ddc28e61919dfcb6152..18311bde1e1bfddb5b1ee566e75da327f79c5f9e 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 1a8e24d07213544a52c16c9ebd89c98cc9b65db9..9cf36e28a2c071bbdab359087074b3ebd449fd17 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 077a3241930d68cfb2825384f7c7ee9400cbd23a..da81a72076fb37774a80ec9ec4b8acbd06971e8a 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/box/localresidual.hh b/dumux/implicit/box/localresidual.hh index 3dcb322ecb997d5af2fa2e26630eddd245b85b6f..bac97a199bbc6f7d9182bfff7ce2044558d98700 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/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh index 606a9e503d609359a6c08b35b4b5c9063a6f1ed2..ef3a5f2a91689a549c8725abf974d012fc27aa53 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 c09efacc7cfa010884b16827877f1629683f2bdc..719dc8a9e5661439cc809d8d212cf8bcdfe9e48e 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.dirichletWrapper(element, scvf); ghostVolVars_[fIdx].update(dirichletValues, problem, diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh index 05d73ea4d9b25f43fb73770b14e58f4f3aad8d93..3b2b386e953ce0911e9d5ea004bc0f38510018a3 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); @@ -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/cellcentered/mpfa/localresidual.hh b/dumux/implicit/cellcentered/mpfa/localresidual.hh index 6fd5394927b90517dfd8277a501251b19f94edd4..a91ea4647079d392868af7d9f9b7216ae2f67a63 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 383200e7e727466d22c795c0179cdd1450d61e4c..c3508d1e38d07b1da7ef44b8f6e5a74b08bffee8 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.dirichletWrapper(element, scvf); (*this)[indexInVariables].update(dirichletValues, problem, diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh index ee68de241760a5110e25b0d41dc30a6a9a6676eb..b62f59863452cfa9f133ea1e6f224bd58c3750fc 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/model.hh b/dumux/implicit/model.hh index fb14955585f92ff413c873e034b84387eca4f399..141d2a76798003993ad399f0850093adfb5e09c3 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 5bbeb00281b05dda994b3f68c7a19cbc6cb53196..89cab6bb86e54e1c347adfc65f294349e54df118 100644 --- a/dumux/implicit/problem.hh +++ b/dumux/implicit/problem.hh @@ -26,12 +26,280 @@ #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().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(...); + + typedef decltype(check(0)) type; + +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; +}; + +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 @@ -150,106 +418,172 @@ 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 values The boundary types for the conservation equations - * \param vertex The vertex for which the boundary type is set + * \param element The element + * \param scv The sub control volume + * + * \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 values The boundary types for the conservation equations - * \param globalPos The position of the finite volume in global coordinates + * \param element The element + * \param scvf The sub control volume face + * + * \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 + { + 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); + break; + } + 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 { - // Throw an exception (there is no reasonable default value - // for Dirichlet conditions) - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide " - "a boundaryTypes() method."); + BoundaryTypes bcTypes; + asImp_().boundaryTypesAtPos(bcTypes, scvf.ipGlobal()); + 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 + { + 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 { - // 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."); + 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; } /*! @@ -260,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 + { + 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 { - // 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."); + PrimaryVariables values(0); + asImp_().neumannAtPos(values, scvf.ipGlobal()); + return values; } /*! @@ -312,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 @@ -324,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; } /*! @@ -417,65 +809,73 @@ 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()); } + 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 + { + FVElementGeometry fvGeometry(model_.globalFvGeometry()); + fvGeometry.bind(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; + } /*! - * \brief Evaluate the initial value for a control volume. + * \brief Returns the initial phase state for a sub 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) + * \param element The element + * \param scv The sub control volume * - * For this method, the \a values parameter stores primary variables. + * \todo replace by initialPhasePresence after release of next */ - PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, int> + initialPhasePresenceWrapper(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."); + return asImp_().initialPhasePresence(element, scv); } - - /*! - * \brief Evaluate the initial phase state inside a control volume. - * - * \param vertex The vertex - * \param vIdxGlobal The global index of the vertex - * \param globalPos The global position - */ - int initialPhasePresence(const Vertex &vertex, - int &vIdxGlobal, - const GlobalPosition &globalPos) const + template + 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> + initialPhasePresenceWrapper(const Element& element, const SubControlVolume& scv) const { - // forward to generic interface - return asImp_().initialPhasePresenceAtPos(globalPos); + int dofIndex = scv.dofIndex(); + return asImp_().initialPhasePresence(element.template subEntity(scv.index()), dofIndex, scv.dofPosition()); } - - /*! - * \brief Evaluate the initial value for a control volume. - * - * \param globalPos The global position - */ - void initialPhasePresenceAtPos(const GlobalPosition &globalPos) const + template + typename std::enable_if_t::value, int> + initialPhasePresenceWrapper(const Element& element, const SubControlVolume& scv) const { - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide a initialPhasePresenceAtPos() method."); + return asImp_().initialPhasePresenceAtPos(scv.dofPosition()); } /*! diff --git a/dumux/porousmediumflow/2p2c/implicit/model.hh b/dumux/porousmediumflow/2p2c/implicit/model.hh index 02c151390e78d097d2882f7ff8061bdd4d39a981..b3162ed8ecb4c2c0124478604fbde4da10c3950c 100644 --- a/dumux/porousmediumflow/2p2c/implicit/model.hh +++ b/dumux/porousmediumflow/2p2c/implicit/model.hh @@ -188,12 +188,12 @@ 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(); 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 650a5f10fcda9d242a0f023e14038fff8925bf3c..7fa8131f0e843c68ef79105bfa20e063fd7d2712 100644 --- a/dumux/porousmediumflow/2pncmin/implicit/model.hh +++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh @@ -195,12 +195,12 @@ 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(); 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/co2/implicit/volumevariables.hh b/dumux/porousmediumflow/co2/implicit/volumevariables.hh index ff328761465fd08f9a372e32b9312981cc956316..a217959c2f4e45cd33060f4ecf89355b8a40ebb7 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 = BaseClassType::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 e8b9a87d2f140c15749c07fb63a47cee66cb4723..7bcd93c290df519eeaf65b75bb443f7577325b64 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.initialPhasePresenceWrapper(element, scv); } } diff --git a/dumux/porousmediumflow/implicit/velocityoutput.hh b/dumux/porousmediumflow/implicit/velocityoutput.hh index 8f386f10425bf0e767c96c664a50cd7360ead310..71c8ed97281f60f0ec9cc03772b14bdfe504d9e1 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/dumux/porousmediumflow/mpnc/implicit/localresidual.hh b/dumux/porousmediumflow/mpnc/implicit/localresidual.hh index 3e7fd6e8cac0a6d1dcedf4985683bf084b9565cd..4ff2eb94b0d4a0a127a1da035a5561c709ee6fbb 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/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh index 2335de51510212edfb838e4f4f2377dbdb8f17f3..7300238467d85ea37675df307bcef98faa9207a4 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; } // \} diff --git a/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh b/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh index 155955e493700cec3bef0e14d4d9c7f68e481f4d..41233e9c5176f8fb78ec46e0a1e9467b4191cedd 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 diff --git a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh index 5f074bd09d47b3d0886c0fd35995743aee0f06db..844878803523bcb3a82d0e3219b98751b3bee09e 100644 --- a/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh +++ b/test/porousmediumflow/co2/implicit/heterogeneousspatialparameters.hh @@ -73,15 +73,18 @@ 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; 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 }; 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: @@ -93,8 +96,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: @@ -150,12 +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 FVElementGeometry &fvGeometry, - int scvIdx) 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_; @@ -172,9 +174,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 +196,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_; }