26 #ifndef EWOMS_VCFV_GRADIENT_CALCULATOR_HH
27 #define EWOMS_VCFV_GRADIENT_CALCULATOR_HH
33 #include <dune/common/fvector.hh>
34 #include <dune/common/version.hh>
36 #include <dune/geometry/type.hh>
37 #include <dune/localfunctions/lagrange/pqkfactory.hh>
38 #include <dune/common/fvector.hh>
50 template<
class TypeTag>
54 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
55 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
56 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
58 enum { dim = GridView::dimension };
63 enum { maxDof = (2 << dim) };
64 enum { maxFap = maxDof };
66 typedef typename GridView::ctype CoordScalar;
67 typedef Dune::FieldVector<Scalar, dim> DimVector;
69 typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
70 typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
71 typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
72 typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
81 template <
bool prepareValues = true,
bool prepareGradients = true>
82 void prepare(
const ElementContext &elemCtx,
int timeIdx)
85 ParentType::template prepare<prepareValues, prepareGradients>(elemCtx, timeIdx);
89 const auto &stencil = elemCtx.stencil(timeIdx);
91 const LocalFiniteElement& localFE = feCache_.get(elemCtx.element().type());
92 localFiniteElement_ = &localFE;
95 for (
int faceIdx = 0; faceIdx < stencil.numInteriorFaces(); ++faceIdx) {
96 const auto &localFacePos = stencil.interiorFace(faceIdx).localPos();
101 localFE.localBasis().evaluateFunction(localFacePos, p1Value_[faceIdx]);
103 if (prepareGradients) {
105 std::vector<ShapeJacobian> localGradient;
106 localFE.localBasis().evaluateJacobian(localFacePos, localGradient);
111 const auto &geom = elemCtx.element().geometry();
112 const auto &jacInvT =
113 geom.jacobianInverseTransposed(localFacePos);
115 int numVertices = elemCtx.numDof(timeIdx);
116 for (
int vertIdx = 0; vertIdx < numVertices; vertIdx++) {
117 jacInvT.mv(localGradient[vertIdx][0],
118 p1Gradient_[faceIdx][vertIdx]);
136 template <
class QuantityCallback,
class Dummy=
int>
138 typename std::enable_if<!
GET_PROP_VALUE(TypeTag, UseTwoPointGradients),
140 const QuantityCallback& quantityCallback)
const
141 ->
typename std::remove_reference<typename QuantityCallback::ResultType>::type
143 static_assert(std::is_integral<Dummy>::value,
144 "The 'Dummy' template parameter must _not_ be specified explicitly."
145 "It is only required to conditionally disable this method!");
147 typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
148 typedef typename std::remove_const<QuantityConstType>::type QuantityType;
152 QuantityType value(0.0);
154 for (
int vertIdx = 0; vertIdx < elemCtx.numDof(0); ++vertIdx) {
155 tmp = quantityCallback(vertIdx);
156 tmp *= p1Value_[fapIdx][vertIdx];
162 template <
class QuantityCallback,
class Dummy =
int>
164 typename std::enable_if<
GET_PROP_VALUE(TypeTag, UseTwoPointGradients),
166 const QuantityCallback& quantityCallback)
const
167 -> decltype(ParentType::calculateValue(elemCtx, fapIdx, quantityCallback))
169 static_assert(std::is_integral<Dummy>::value,
170 "The 'Dummy' template parameter must _not_ be specified explicitly."
171 "It is only required to conditionally disable this method!");
173 return ParentType::calculateValue(elemCtx, fapIdx, quantityCallback);
187 template <
class QuantityCallback,
class EvalDimVector,
class Dummy =
int>
189 const ElementContext &elemCtx,
190 typename std::enable_if<!
GET_PROP_VALUE(TypeTag, UseTwoPointGradients),
192 const QuantityCallback& quantityCallback)
const
194 static_assert(std::is_integral<Dummy>::value,
195 "The 'Dummy' template parameter must _not_ be specified explicitly."
196 "It is only required to conditionally disable this method!");
201 for (
int vertIdx = 0; vertIdx < elemCtx.numDof(0); ++vertIdx) {
202 Scalar dofVal = quantityCallback(vertIdx);
204 auto tmp = p1Gradient_[fapIdx][vertIdx];
210 template <
class QuantityCallback,
class EvalDimVector,
class Dummy=
int>
212 const ElementContext &elemCtx,
213 typename std::enable_if<
GET_PROP_VALUE(TypeTag, UseTwoPointGradients),
215 const QuantityCallback& quantityCallback)
const
217 static_assert(std::is_integral<Dummy>::value,
218 "The 'Dummy' template parameter must _not_ be specified explicitly."
219 "It is only required to conditionally disable this method!");
221 return ParentType::calculateGradient(quantityGrad, elemCtx, fapIdx, quantityCallback);
238 template <
class QuantityCallback>
241 const QuantityCallback &quantityCallback)
242 -> decltype(ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback))
243 {
return ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback); }
259 template <
class QuantityCallback,
class EvalDimVector>
261 const ElementContext &elemCtx,
263 const QuantityCallback &quantityCallback)
const
264 { ParentType::calculateBoundaryGradient(quantityGrad, elemCtx, fapIdx, quantityCallback); }
270 static LocalFiniteElementCache feCache_;
272 const LocalFiniteElement* localFiniteElement_;
273 std::vector<Dune::FieldVector<Scalar, 1>> p1Value_[maxFap];
274 DimVector p1Gradient_[maxFap][maxDof];
277 template<
class TypeTag>
278 typename VcfvGradientCalculator<TypeTag>::LocalFiniteElementCache
279 VcfvGradientCalculator<TypeTag>::feCache_;
static LocalFiniteElementCache & localFiniteElementCache()
Definition: vcfvgradientcalculator.hh:266
auto calculateValue(const ElementContext &elemCtx, typename std::enable_if< GET_PROP_VALUE(TypeTag, UseTwoPointGradients), Dummy >::type fapIdx, const QuantityCallback &quantityCallback) const -> decltype(ParentType::calculateValue(elemCtx, fapIdx, quantityCallback))
Definition: vcfvgradientcalculator.hh:163
void calculateGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, typename std::enable_if<!GET_PROP_VALUE(TypeTag, UseTwoPointGradients), Dummy >::type fapIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point. ...
Definition: vcfvgradientcalculator.hh:188
void calculateBoundaryGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, int fapIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point on the boundary...
Definition: vcfvgradientcalculator.hh:260
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
void prepare(const ElementContext &elemCtx, int timeIdx)
Precomputes the common values to calculate gradients and values of quantities at any flux approximati...
Definition: vcfvgradientcalculator.hh:82
auto calculateValue(const ElementContext &elemCtx, typename std::enable_if<!GET_PROP_VALUE(TypeTag, UseTwoPointGradients), Dummy >::type fapIdx, const QuantityCallback &quantityCallback) const -> typename std::remove_reference< typename QuantityCallback::ResultType >::type
Calculates the value of an arbitrary quantity at any interior flux approximation point.
Definition: vcfvgradientcalculator.hh:137
Definition: baseauxiliarymodule.hh:35
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
This class calculates gradients of arbitrary quantities at flux integration points for the vertex cen...
Definition: vcfvgradientcalculator.hh:51
auto calculateBoundaryValue(const ElementContext &elemCtx, int fapIdx, const QuantityCallback &quantityCallback) -> decltype(ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback))
Calculates the value of an arbitrary quantity at any flux approximation point on the grid boundary...
Definition: vcfvgradientcalculator.hh:239
void calculateGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, typename std::enable_if< GET_PROP_VALUE(TypeTag, UseTwoPointGradients), Dummy >::type fapIdx, const QuantityCallback &quantityCallback) const
Definition: vcfvgradientcalculator.hh:211
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:44
Declares the basic properties used by the common infrastructure of the vertex-centered finite volume ...