25#ifndef DUMUX_PNM_UTILITES_HH
26#define DUMUX_PNM_UTILITES_HH
38template<
class Gr
idVariables,
class SolutionVector>
41 using Scalar =
typename GridVariables::VolumeVariables::PrimaryVariables::value_type;
42 using VolumeVariables =
typename GridVariables::VolumeVariables;
44 struct AveragedQuantityInfo
46 std::function<Scalar(
const VolumeVariables&)> quantity;
47 std::function<Scalar(
const VolumeVariables&)> weight;
54 const SolutionVector& sol)
55 : gridVariables_(gridVariables)
60 std::function<Scalar(
const VolumeVariables&)>&& w,
61 const std::string& name)
63 averagedQuantityInfo_.push_back(AveragedQuantityInfo{q, w, name});
64 averagedQuantity_.push_back(Scalar());
67 void eval(
const std::vector<std::size_t>& dofsToNeglect = std::vector<std::size_t>())
69 for (
auto& q : averagedQuantity_)
72 auto fvGeometry =
localView(problem_().gridGeometry());
73 auto elemVolVars =
localView(gridVariables_.curGridVolVars());
74 std::vector<bool> poreVisited(problem_().gridGeometry().numDofs(),
false);
75 std::vector<Scalar> weights(averagedQuantityInfo_.size(), 0.0);
77 for (
const auto& element : elements(problem_().gridGeometry().gridView()))
79 fvGeometry.bind(element);
80 elemVolVars.bind(element, fvGeometry, sol_);
82 for (
int scvIdx = 0; scvIdx < fvGeometry.numScv(); ++scvIdx)
84 static constexpr auto dofCodim = std::decay_t<
decltype(problem_().gridGeometry().gridView())>::dimension;
85 const int dofIdxGlobal = problem_().gridGeometry().vertexMapper().subIndex(element, scvIdx, dofCodim);
87 if (poreVisited[dofIdxGlobal])
89 else if (!dofsToNeglect.empty() && std::any_of(dofsToNeglect.begin(), dofsToNeglect.end(), [&](
int dofIdx){ return dofIdx == dofIdxGlobal; }))
93 const auto& volVars = elemVolVars[scvIdx];
94 for (
int i = 0; i < averagedQuantityInfo_.size(); ++i)
96 const Scalar weight = averagedQuantityInfo_[i].weight(volVars);
97 averagedQuantity_[i] += averagedQuantityInfo_[i].quantity(volVars) * weight;
100 poreVisited[dofIdxGlobal] =
true;
105 for (
int i = 0; i < averagedQuantityInfo_.size(); ++i)
106 averagedQuantity_[i] /= weights[i];
111 for (
int i = 0; i < averagedQuantityInfo_.size(); ++i)
112 if (averagedQuantityInfo_[i].name == name)
113 return averagedQuantity_[i];
115 DUNE_THROW(Dune::InvalidStateException, name +
" not found");
120 std::vector<AveragedQuantityInfo> averagedQuantityInfo_;
121 std::vector<Scalar> averagedQuantity_;
123 const auto& problem_()
124 {
return gridVariables_.curGridVolVars().problem(); }
126 const GridVariables& gridVariables_;
127 const SolutionVector& sol_;
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Definition: discretization/porenetwork/fvelementgeometry.hh:33
Calculates averaged values of the network solution.
Definition: utilities.hh:40
AveragedValues(const GridVariables &gridVariables, const SolutionVector &sol)
Definition: utilities.hh:53
void addAveragedQuantity(std::function< Scalar(const VolumeVariables &)> &&q, std::function< Scalar(const VolumeVariables &)> &&w, const std::string &name)
Definition: utilities.hh:59
void eval(const std::vector< std::size_t > &dofsToNeglect=std::vector< std::size_t >())
Definition: utilities.hh:67
const Scalar & operator[](const std::string &name) const
Definition: utilities.hh:109