Skip to content
GitLab
Commit 2d684fe1 authored by Martin Schneider's avatar Martin Schneider
Browse files

[md][stokesdarcy][box] Introduce separate class to calculate projection

parent 20265b63
Hide whitespace changes
Inline Side-by-side
dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh
View file @ 2d684fe1
......@@ -34,6 +34,7 @@
#include <dumux/freeflow/navierstokes/staggered/velocitygradients.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include "projection.hh"
#include <optional>
namespace Dumux {
......@@ -73,6 +74,9 @@ class StokesDarcyCouplingDataBoxBase : public StokesDarcyCouplingDataImplementat
using ForchheimersLaw = ForchheimersLawImplementation<SubDomainTypeTag<porousMediumIdx>, GridGeometry<porousMediumIdx>::discMethod>;
using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
using Projector = Projection<MDTraits, CouplingManager>;
public:
StokesDarcyCouplingDataBoxBase(const CouplingManager& couplingmanager): ParentType(couplingmanager) {}
......@@ -100,7 +104,7 @@ public:
auto pressure = [darcyPhaseIdx](const auto& elemVolVars, const auto& scv)
{ return elemVolVars[scv].pressure(darcyPhaseIdx); };
Scalar momentumFlux = this->couplingManager().calculateProjection(element, scvf, pressure);
Scalar momentumFlux = Projector::calculateProjection(this->couplingManager(), element, scvf, pressure);
// normalize pressure
if(getPropValue<SubDomainTypeTag<freeFlowIdx>, Properties::NormalizePressure>())
......@@ -190,7 +194,7 @@ public:
VelocityVector velocity(0.0); // velocity darcy // density darcy
Scalar intersectionLength = 0.0; // (total)intersection length, could differ from scvf length
const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf);
const auto& stokesContext = this->couplingManager().stokesCouplingContext();
static const bool enableGravity = getParamFromGroup<bool>(this->couplingManager().problem(porousMediumIdx).paramGroup(), "Problem.EnableGravity");
// iteration over the different coupling facets
......@@ -295,6 +299,8 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
using Projector = Projection<MDTraits, CouplingManager>;
public:
using ParentType::ParentType;
using ParentType::couplingPhaseIdx;
......@@ -309,7 +315,7 @@ public:
const SubControlVolumeFace<porousMediumIdx>& scvf) const
{
const auto darcyPhaseIdx = couplingPhaseIdx(porousMediumIdx);
const auto& darcyContext = this->couplingManager().darcyCouplingContextVector(element, scvf);
const auto& darcyContext = this->couplingManager().darcyCouplingContext();
Scalar flux = 0.0;
for(const auto& data : darcyContext)
......@@ -341,7 +347,7 @@ public:
const ElementFaceVariables<freeFlowIdx>& stokesElemFaceVars,
const SubControlVolumeFace<freeFlowIdx>& scvf) const
{
const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf);
const auto& stokesContext = this->couplingManager().stokesCouplingContext();
Scalar flux = 0.0;
const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf() * scvf.directionSign();
......@@ -378,7 +384,7 @@ public:
const SubControlVolumeFace<porousMediumIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& darcyContext = this->couplingManager().darcyCouplingContextVector(element, scvf);
const auto& darcyContext = this->couplingManager().darcyCouplingContext();
Scalar flux = 0.0;
for(const auto& data : darcyContext)
......@@ -396,7 +402,7 @@ public:
//Calculate the projected temperature value for the stokes face
auto temp = [](const auto& elemVolVars, const auto& scv)
{ return elemVolVars[scv].temperature(); };
Scalar interfaceTemperature = this->couplingManager().calculateProjection(stokesScvf, element, darcyElemVolVars, temp);
Scalar interfaceTemperature = Projector::calculateProjection(this->couplingManager(), stokesScvf, element, darcyElemVolVars, temp);
const bool insideIsUpstream = velocity > 0.0;
flux += -1*energyFlux_(data.fvGeometry,
......@@ -424,13 +430,13 @@ public:
const SubControlVolumeFace<freeFlowIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf);
const auto& stokesContext = this->couplingManager().stokesCouplingContext();
//Calculate the projected temperature value for the stokes face
auto temp = [](const auto& elemVolVars, const auto& scv)
{ return elemVolVars[scv].temperature(); };
Scalar interfaceTemperature = this->couplingManager().calculateProjection(element, scvf, temp);
Scalar interfaceTemperature = Projector::calculateProjection(this->couplingManager(), element, scvf, temp);
Scalar flux = 0.0;
for (const auto& data : stokesContext)
......@@ -562,6 +568,8 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
using MolecularDiffusionType = GetPropType<SubDomainTypeTag<freeFlowIdx>, Properties::MolecularDiffusionType>;
using Projector = Projection<MDTraits, CouplingManager>;
public:
using ParentType::ParentType;
using ParentType::couplingPhaseIdx;
......@@ -577,7 +585,7 @@ public:
const SubControlVolumeFace<porousMediumIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& darcyContext = this->couplingManager().darcyCouplingContextVector(element, scvf);
const auto& darcyContext = this->couplingManager().darcyCouplingContext();
NumEqVector flux(0.0);
for(const auto& data : darcyContext)
......@@ -599,7 +607,7 @@ public:
auto value = [&compIdx](const auto& elemVolVars, const auto& scv)
{ return massOrMoleFraction(elemVolVars[scv], referenceSystemFormulation, couplingPhaseIdx(porousMediumIdx), compIdx); };
return this->couplingManager().calculateProjection(stokesScvf, element, darcyElemVolVars, value);
return Projector::calculateProjection(this->couplingManager(), stokesScvf, element, darcyElemVolVars, value);
};
// Division by scvf.area() is needed, because the final flux results from multiplication with scvf.area()
......@@ -631,7 +639,7 @@ public:
const SubControlVolumeFace<freeFlowIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf);
const auto& stokesContext = this->couplingManager().stokesCouplingContext();
//Calculate the projected massOrMoleFraction value for the stokes face
auto interfaceMassOrMoleFraction = [this, &element, &scvf](int compIdx)
......@@ -639,7 +647,7 @@ public:
auto value = [&compIdx](const auto& elemVolVars, const auto& scv)
{ return massOrMoleFraction(elemVolVars[scv], referenceSystemFormulation, couplingPhaseIdx(porousMediumIdx), compIdx); };
return this->couplingManager().calculateProjection(element, scvf, value);
return Projector::calculateProjection(this->couplingManager(), element, scvf, value);
};
NumEqVector flux(0.0);
......@@ -684,7 +692,7 @@ public:
const SubControlVolumeFace<porousMediumIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& darcyContext = this->couplingManager().darcyCouplingContextVector(element, scvf);
const auto& darcyContext = this->couplingManager().darcyCouplingContext();
Scalar flux = 0.0;
for(const auto& data : darcyContext)
......@@ -706,13 +714,13 @@ public:
auto value = [&compIdx](const auto& elemVolVars, const auto& scv)
{ return massOrMoleFraction(elemVolVars[scv], referenceSystemFormulation, couplingPhaseIdx(porousMediumIdx), compIdx); };
return this->couplingManager().calculateProjection(stokesScvf, element, darcyElemVolVars, value);
return Projector::calculateProjection(this->couplingManager(), stokesScvf, element, darcyElemVolVars, value);
};
//Calculate the projected temperature value for the stokes face
auto temp = [](const auto& elemVolVars, const auto& scv)
{ return elemVolVars[scv].temperature(); };
Scalar interfaceTemperature = this->couplingManager().calculateProjection(stokesScvf, element, darcyElemVolVars, temp);
Scalar interfaceTemperature = Projector::calculateProjection(this->couplingManager(), stokesScvf, element, darcyElemVolVars, temp);
// Division by scvf.area() is needed, because the final flux results from multiplication with scvf.area()
flux += -1*energyFlux_(porousMediumIdx,
......@@ -745,7 +753,7 @@ public:
const SubControlVolumeFace<freeFlowIdx>& scvf,
const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
{
const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf);
const auto& stokesContext = this->couplingManager().stokesCouplingContext();
//Calculate the projected massOrMoleFraction value for the stokes face
auto interfaceMassOrMoleFraction = [this, &element, &scvf](int compIdx)
......@@ -753,14 +761,14 @@ public:
auto value = [&compIdx](const auto& elemVolVars, const auto& scv)
{ return massOrMoleFraction(elemVolVars[scv], referenceSystemFormulation, couplingPhaseIdx(porousMediumIdx), compIdx); };
return this->couplingManager().calculateProjection(element, scvf, value);
return Projector::calculateProjection(this->couplingManager(), element, scvf, value);
};
//Calculate the projected temperature value for the stokes face
auto temp = [](const auto& elemVolVars, const auto& scv)
{ return elemVolVars[scv].temperature(); };
Scalar interfaceTemperature = this->couplingManager().calculateProjection(element, scvf, temp);
Scalar interfaceTemperature = Projector::calculateProjection(this->couplingManager(), element, scvf, temp);
Scalar flux = 0.0;
const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf() * scvf.directionSign();
......
dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh
View file @ 2d684fe1
......@@ -43,28 +43,6 @@
namespace Dumux {
namespace Detail {
enum class ProjectionMethod
{
L2Projection, AreaWeightedDofEvaluation
};
template <typename T>
using ProjectionMethodDetector = decltype(std::declval<T>().projectionMethod());
template<class T>
static constexpr bool hasProjectionMethod()
{ return Dune::Std::is_detected<ProjectionMethodDetector, T>::value; }
template<class T>
static constexpr ProjectionMethod projectionMethod()
{
if constexpr (hasProjectionMethod<T>())
return T::projectionMethod();
else
return ProjectionMethod::L2Projection;
}
// Each context object contains the data related to one coupling facet
template <class MDTraits, typename CouplingFacetGeometry>
struct StokesCouplingContext
......@@ -144,8 +122,6 @@ public:
static constexpr auto freeFlowIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::freeFlowIdx;
static constexpr auto porousMediumIdx = FreeFlowPorousMediumCouplingManagerBase<MDTraits>::porousMediumIdx;
using ProjectionMethod = Detail::ProjectionMethod;
static constexpr auto projectionMethod = Detail::projectionMethod<MDTraits>();
private:
using SolutionVector = typename MDTraits::SolutionVector;
......@@ -460,10 +436,10 @@ public:
/*!
* \brief Access the coupling context needed for the Darcy domain
*/
const auto& darcyCouplingContextVector(const Element<porousMediumIdx>& element, const SubControlVolumeFace<porousMediumIdx>& scvf) const
const auto& darcyCouplingContext() const
{
if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element))
DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
if (darcyCouplingContext_.empty())
DUNE_THROW(Dune::InvalidStateException, "Coupling context is empty!");
return darcyCouplingContext_;
}
......@@ -471,10 +447,11 @@ public:
/*!
* \brief Access the coupling context needed for the Stokes domain
*/
const auto& stokesCouplingContextVector(const Element<freeFlowIdx>& element, const SubControlVolumeFace<freeFlowIdx>& scvf) const
const auto& stokesCouplingContext() const
{
if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx())
DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
if (stokesCouplingContext_.empty())
DUNE_THROW(Dune::InvalidStateException, "Coupling context is empty!");
return stokesCouplingContext_;
}
......@@ -727,115 +704,14 @@ public:
}
}
// calculate projection of pm solution needed for fluxes of the free-flow residual
template<class Function>
Scalar calculateProjection(const SubControlVolumeFace<freeFlowIdx>& stokesScvf,
const Element<porousMediumIdx>& darcyElement,
const ElementVolumeVariables<porousMediumIdx>& darcyElemVolVars,
Function evalPriVar) const
{
Scalar projection = 0.0;
auto domainI = Dune::index_constant<freeFlowIdx>();
auto fvGeometry = localView(this->problem(porousMediumIdx).gridGeometry());
auto elemVolVars = localView(darcyElemVolVars.gridVolVars());
// integrate darcy pressure over each coupling facet and average
for(const auto& couplingFacet : couplingFacets(domainI, couplingMapper_, stokesScvf.insideScvIdx(), stokesScvf.localFaceIdx()))
{
const auto darcyEIdxI = this->problem(porousMediumIdx).gridGeometry().elementMapper().index(darcyElement);
const auto darcyEIdxJ = couplingFacet.pmEIdx;
const auto& element = this->problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.pmEIdx);
fvGeometry.bind(element);
if(darcyEIdxI == darcyEIdxJ)
{
projection += calculateFacetIntegral(element, fvGeometry, fvGeometry.scvf(couplingFacet.pmScvfIdx), darcyElemVolVars, couplingFacet.geometry, evalPriVar);
}
else
{
elemVolVars.bind(element, fvGeometry, this->curSol()[porousMediumIdx]);
projection += calculateFacetIntegral(element, fvGeometry, fvGeometry.scvf(couplingFacet.pmScvfIdx), elemVolVars, couplingFacet.geometry, evalPriVar);
}
}
projection /= stokesScvf.area();
return projection;
}
// calculate projection of pm solution needed for fluxes of ff residual
template<class Function>
Scalar calculateProjection(const Element<freeFlowIdx>& element,
const SubControlVolumeFace<freeFlowIdx>& scvf,
Function evalPriVar) const
const auto sol() const
{
if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx())
DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center());
Scalar projection = 0.0;
// integrate darcy pressure over each coupling facet and average
for (const auto& data : stokesCouplingContext_)
{
//ToDo Is this if really necessary?
if (scvf.index() == data.stokesScvfIdx)
{
const auto& elemVolVars = *(data.elementVolVars);
const auto& darcyScvf = data.fvGeometry.scvf(data.darcyScvfIdx);
const auto& couplingFacet = couplingMapper_.couplingFacet(data.facetIdx);
projection += calculateFacetIntegral(data.element, data.fvGeometry, darcyScvf, elemVolVars, couplingFacet.geometry, evalPriVar);
}
}
projection /= scvf.area();
return projection;
return this->curSol();
}
template<class CouplingFacetGeometry, class Function>
Scalar calculateFacetIntegral(const Element<porousMediumIdx>& element,
const FVElementGeometry<porousMediumIdx>& fvGeometry,
const SubControlVolumeFace<porousMediumIdx>& scvf,
const ElementVolumeVariables<porousMediumIdx>& elemVolVars,
const CouplingFacetGeometry& facetGeometry,
Function evalPriVar) const
const auto& couplingMapper() const
{
Scalar facetProjection = 0.0;
if constexpr (projectionMethod == ProjectionMethod::L2Projection)
{
const auto& localBasis = fvGeometry.feLocalBasis();
// do second order integration as box provides linear functions
static constexpr int darcyDim = GridGeometry<porousMediumIdx>::GridView::dimension;
const auto& rule = Dune::QuadratureRules<Scalar, darcyDim-1>::rule(facetGeometry.type(), 2);
for (const auto& qp : rule)
{
const auto& ipLocal = qp.position();
const auto& ipGlobal = facetGeometry.global(ipLocal);
const auto& ipElementLocal = element.geometry().local(ipGlobal);
std::vector<Dune::FieldVector<Scalar, 1>> shapeValues;
localBasis.evaluateFunction(ipElementLocal, shapeValues);
Scalar value = 0.0;
for (const auto& scv : scvs(fvGeometry))
value += evalPriVar(elemVolVars, scv)*shapeValues[scv.indexInElement()][0];
facetProjection += value*facetGeometry.integrationElement(qp.position())*qp.weight();
}
}
else if constexpr (projectionMethod == ProjectionMethod::AreaWeightedDofEvaluation)
{
facetProjection = facetGeometry.volume()*evalPriVar(elemVolVars, fvGeometry.scv(scvf.insideScvIdx()));
}
else
{
DUNE_THROW(Dune::NotImplemented, "Unkown projection method!");
}
return facetProjection;
return couplingMapper_;
}
const auto& couplingFacet(std::size_t idx) const
......
dumux/multidomain/boundary/stokesdarcy/box/projection.hh 0 → 100644
View file @ 2d684fe1
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup StokesDarcyCoupling
*/
#ifndef DUMUX_STOKES_DARCY_PROJECTION_HH
#define DUMUX_STOKES_DARCY_PROJECTION_HH
#include <type_traits>
#include <dumux/common/properties.hh>
#include "couplingmanager.hh"
namespace Dumux {
namespace Detail {
enum class ProjectionMethod
{
L2Projection, AreaWeightedDofEvaluation
};
template <typename T>
using ProjectionMethodDetector = decltype(std::declval<T>().projectionMethod());
template<class T>
static constexpr bool hasProjectionMethod()
{ return Dune::Std::is_detected<ProjectionMethodDetector, T>::value; }
template<class T>
static constexpr ProjectionMethod projectionMethod()
{
if constexpr (hasProjectionMethod<T>())
return T::projectionMethod();
else
return ProjectionMethod::L2Projection;
}
}
template<class MDTraits, class CouplingManager>
class Projection
{
using Scalar = typename MDTraits::Scalar;
static constexpr auto freeFlowIdx = CouplingManager::freeFlowIdx;
static constexpr auto porousMediumIdx = CouplingManager::porousMediumIdx;
// the sub domain type tags
template<std::size_t id>
using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
using ProjectionMethod = Detail::ProjectionMethod;
static constexpr auto projectionMethod = Detail::projectionMethod<MDTraits>();
public:
// calculate projection of pm solution needed for fluxes of the free-flow residual
template<class Function>
static Scalar calculateProjection(const CouplingManager& couplingManager,
const SubControlVolumeFace<freeFlowIdx>& stokesScvf,
const Element<porousMediumIdx>& darcyElement,
const ElementVolumeVariables<porousMediumIdx>& darcyElemVolVars,
Function evalPriVar)
{
Scalar projection = 0.0;
auto domainI = Dune::index_constant<freeFlowIdx>();
auto fvGeometry = localView(couplingManager.problem(porousMediumIdx).gridGeometry());
auto elemVolVars = localView(darcyElemVolVars.gridVolVars());
// integrate darcy pressure over each coupling facet and average
for(const auto& couplingFacet : couplingFacets(domainI, couplingManager.couplingMapper(), stokesScvf.insideScvIdx(), stokesScvf.localFaceIdx()))
{
const auto darcyEIdxI = couplingManager.problem(porousMediumIdx).gridGeometry().elementMapper().index(darcyElement);
const auto darcyEIdxJ = couplingFacet.pmEIdx;
const auto& element = couplingManager.problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingFacet.pmEIdx);
fvGeometry.bind(element);
if(darcyEIdxI == darcyEIdxJ)
{
projection += calculateFacetIntegral(element, fvGeometry, fvGeometry.scvf(couplingFacet.pmScvfIdx), darcyElemVolVars, couplingFacet.geometry, evalPriVar);
}
else
{
elemVolVars.bind(element, fvGeometry, couplingManager.sol()[porousMediumIdx]);
projection += calculateFacetIntegral(element, fvGeometry, fvGeometry.scvf(couplingFacet.pmScvfIdx), elemVolVars, couplingFacet.geometry, evalPriVar);
}
}
projection /= stokesScvf.area();
return projection;
}
// calculate projection of pm solution needed for fluxes of ff residual
template<class Function>
static Scalar calculateProjection(const CouplingManager& couplingManager,
const Element<freeFlowIdx>& element,
const SubControlVolumeFace<freeFlowIdx>& scvf,
Function evalPriVar)
{
Scalar projection = 0.0;
// integrate darcy pressure over each coupling facet and average
for (const auto& data : couplingManager.stokesCouplingContext())
{
//ToDo Is this if really necessary?
if (scvf.index() == data.stokesScvfIdx)
{
const auto& elemVolVars = *(data.elementVolVars);
const auto& darcyScvf = data.fvGeometry.scvf(data.darcyScvfIdx);
const auto& couplingFacet = couplingManager.couplingMapper().couplingFacet(data.facetIdx);
projection += calculateFacetIntegral(data.element, data.fvGeometry, darcyScvf, elemVolVars, couplingFacet.geometry, evalPriVar);
}
}
projection /= scvf.area();
return projection;
}
template<class CouplingFacetGeometry, class Function>
static Scalar calculateFacetIntegral(const Element<porousMediumIdx>& element,
const FVElementGeometry<porousMediumIdx>& fvGeometry,
const SubControlVolumeFace<porousMediumIdx>& scvf,
const ElementVolumeVariables<porousMediumIdx>& elemVolVars,
const CouplingFacetGeometry& facetGeometry,
Function evalPriVar)
{
Scalar facetProjection = 0.0;
if constexpr (projectionMethod == ProjectionMethod::L2Projection)
{
const auto& localBasis = fvGeometry.feLocalBasis();
// do second order integration as box provides linear functions
static constexpr int darcyDim = GridGeometry<porousMediumIdx>::GridView::dimension;
const auto& rule = Dune::QuadratureRules<Scalar, darcyDim-1>::rule(facetGeometry.type(), 2);
for (const auto& qp : rule)
{
const auto& ipLocal = qp.position();
const auto& ipGlobal = facetGeometry.global(ipLocal);
const auto& ipElementLocal = element.geometry().local(ipGlobal);
std::vector<Dune::FieldVector<Scalar, 1>> shapeValues;
localBasis.evaluateFunction(ipElementLocal, shapeValues);
Scalar value = 0.0;
for (const auto& scv : scvs(fvGeometry))
value += evalPriVar(elemVolVars, scv)*shapeValues[scv.indexInElement()][0];
facetProjection += value*facetGeometry.integrationElement(qp.position())*qp.weight();
}
}
else if constexpr (projectionMethod == ProjectionMethod::AreaWeightedDofEvaluation)
{
facetProjection = facetGeometry.volume()*evalPriVar(elemVolVars, fvGeometry.scv(scvf.insideScvIdx()));
}
else
{
DUNE_THROW(Dune::NotImplemented, "Unkown projection method!");
}
return facetProjection;
}
};
} // end namespace Dumux
#endif // DUMUX_STOKES_DARCY_PROJECTION_HH
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment