28 #ifndef EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH 29 #define EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH 33 #include <opm/common/Unused.hpp> 35 #include <dune/common/fvector.hh> 38 template<
class TypeTag>
47 template<
class TypeTag>
50 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
51 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
52 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
53 typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
54 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
56 enum { dim = GridView::dimension };
57 enum { dimWorld = GridView::dimensionworld };
61 enum { maxFap = 2 << dim };
63 typedef Opm::MathToolbox<Evaluation> Toolbox;
64 typedef Dune::FieldVector<Scalar, dim> DimVector;
65 typedef Dune::FieldVector<Evaluation, dim> EvalDimVector;
82 template <
bool prepareValues = true,
bool prepareGradients = true>
83 void prepare(
const ElementContext& elemCtx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED)
97 template <
class QuantityCallback>
100 const QuantityCallback& quantityCallback)
const 101 ->
typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
103 typedef decltype(quantityCallback.operator()(0)) RawReturnType;
104 typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
106 typedef Opm::MathToolbox<ReturnType> Toolbox;
108 Scalar interiorDistance;
109 Scalar exteriorDistance;
110 computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
112 const auto& face = elemCtx.stencil(0).interiorFace(fapIdx);
113 auto i = face.interiorIndex();
114 auto j = face.exteriorIndex();
115 auto focusDofIdx = elemCtx.focusDofIndex();
119 if (i == focusDofIdx)
120 value = quantityCallback(i)*interiorDistance;
122 value = Toolbox::value(quantityCallback(i))*interiorDistance;
124 if (j == focusDofIdx)
125 value += quantityCallback(j)*exteriorDistance;
127 value += Toolbox::value(quantityCallback(j))*exteriorDistance;
129 value /= interiorDistance + exteriorDistance;
145 template <
class QuantityCallback>
148 const QuantityCallback& quantityCallback)
const 149 ->
typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
151 typedef decltype(quantityCallback.operator()(0)) RawReturnType;
152 typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
154 typedef decltype(std::declval<ReturnType>()[0]) RawFieldType;
155 typedef typename std::remove_const<typename std::remove_reference<RawFieldType>::type>::type FieldType;
157 typedef Opm::MathToolbox<FieldType> Toolbox;
159 Scalar interiorDistance;
160 Scalar exteriorDistance;
161 computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
163 const auto& face = elemCtx.stencil(0).interiorFace(fapIdx);
164 auto i = face.interiorIndex();
165 auto j = face.exteriorIndex();
166 auto focusDofIdx = elemCtx.focusDofIndex();
170 if (i == focusDofIdx) {
171 value = quantityCallback(i);
172 for (
int k = 0; k < value.size(); ++k)
173 value[k] *= interiorDistance;
176 const auto& dofVal = Toolbox::value(quantityCallback(i));
177 for (
int k = 0; k < dofVal.size(); ++k)
178 value[k] = Toolbox::value(dofVal[k])*interiorDistance;
181 if (j == focusDofIdx) {
182 const auto& dofVal = quantityCallback(j);
183 for (
int k = 0; k < dofVal.size(); ++k)
184 value[k] += dofVal[k]*exteriorDistance;
187 const auto& dofVal = quantityCallback(j);
188 for (
int k = 0; k < dofVal.size(); ++k)
189 value[k] += Toolbox::value(dofVal[k])*exteriorDistance;
192 Scalar totDistance = interiorDistance + exteriorDistance;
193 for (
int k = 0; k < value.size(); ++k)
194 value[k] /= totDistance;
210 template <
class QuantityCallback>
212 const ElementContext& elemCtx,
214 const QuantityCallback& quantityCallback)
const 216 const auto& stencil = elemCtx.stencil(0);
217 const auto& face = stencil.interiorFace(fapIdx);
219 auto i = face.interiorIndex();
220 auto j = face.exteriorIndex();
221 auto focusIdx = elemCtx.focusDofIndex();
223 const auto& interiorPos = stencil.subControlVolume(i).globalPos();
224 const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
229 Toolbox::value(quantityCallback(j))
230 - quantityCallback(i);
232 else if (j == focusIdx) {
235 - Toolbox::value(quantityCallback(i));
239 Toolbox::value(quantityCallback(j))
240 - Toolbox::value(quantityCallback(i));
242 Scalar distSquared = 0.0;
243 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
244 Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
245 distSquared += tmp*tmp;
252 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
253 Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
254 quantityGrad[dimIdx] = deltay*(tmp/distSquared);
272 template <
class QuantityCallback>
274 unsigned fapIdx OPM_UNUSED,
275 const QuantityCallback& quantityCallback)
276 -> decltype(quantityCallback.boundaryValue())
277 {
return quantityCallback.boundaryValue(); }
293 template <
class QuantityCallback>
295 const ElementContext& elemCtx,
297 const QuantityCallback& quantityCallback)
const 299 typedef typename std::decay<decltype(quantityCallback(0))>::type Evaluation;
300 typedef typename std::decay<decltype(quantityCallback.boundaryValue())>::type BoundaryEvaluation;
302 typedef Opm::MathToolbox<Evaluation> Toolbox;
303 typedef Opm::MathToolbox<BoundaryEvaluation> BoundaryToolbox;
305 const auto& stencil = elemCtx.stencil(0);
306 const auto& face = stencil.boundaryFace(faceIdx);
309 if (face.interiorIndex() == elemCtx.focusDofIndex())
310 deltay = quantityCallback.boundaryValue() - quantityCallback(face.interiorIndex());
313 BoundaryToolbox::value(quantityCallback.boundaryValue())
314 - Toolbox::value(quantityCallback(face.interiorIndex()));
316 const auto& boundaryFacePos = face.integrationPos();
317 const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();
319 Scalar distSquared = 0;
320 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
321 Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
322 distSquared += tmp*tmp;
330 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
331 Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
332 quantityGrad[dimIdx] = deltay*(tmp/distSquared);
337 void computeDistances_(Scalar& interiorDistance,
338 Scalar& exteriorDistance,
339 const ElementContext& elemCtx,
340 unsigned fapIdx)
const 342 const auto& stencil = elemCtx.stencil(0);
343 const auto& face = stencil.interiorFace(fapIdx);
347 const auto& normal = face.normal();
348 auto i = face.interiorIndex();
349 auto j = face.exteriorIndex();
350 const auto& interiorPos = stencil.subControlVolume(i).globalPos();
351 const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
352 const auto& integrationPos = face.integrationPos();
354 interiorDistance = 0.0;
355 exteriorDistance = 0.0;
356 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
358 (interiorPos[dimIdx] - integrationPos[dimIdx])
362 (exteriorPos[dimIdx] - integrationPos[dimIdx])
366 interiorDistance = std::sqrt(std::abs(interiorDistance));
367 exteriorDistance = std::sqrt(std::abs(exteriorDistance));
Definition: baseauxiliarymodule.hh:37
void calculateBoundaryGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, unsigned faceIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point on the boundary...
Definition: fvbasegradientcalculator.hh:294
void calculateGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point. ...
Definition: fvbasegradientcalculator.hh:211
static void registerParameters()
Register all run-time parameters for the gradient calculator of the base class of the discretization...
Definition: fvbasegradientcalculator.hh:72
The base class for the element-centered finite-volume discretization scheme.
Definition: fvbasegradientcalculator.hh:39
Declare the properties used by the infrastructure code of the finite volume discretizations.
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:48
void prepare(const ElementContext &elemCtx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Precomputes the common values to calculate gradients and values of quantities at every interior flux ...
Definition: fvbasegradientcalculator.hh:83
auto calculateScalarValue(const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const -> typename std::remove_reference< decltype(quantityCallback.operator()(0))>::type
Calculates the value of an arbitrary scalar quantity at any interior flux approximation point...
Definition: fvbasegradientcalculator.hh:98
auto calculateBoundaryValue(const ElementContext &elemCtx OPM_UNUSED, unsigned fapIdx OPM_UNUSED, const QuantityCallback &quantityCallback) -> decltype(quantityCallback.boundaryValue())
Calculates the value of an arbitrary quantity at any flux approximation point on the grid boundary...
Definition: fvbasegradientcalculator.hh:273
auto calculateVectorValue(const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const -> typename std::remove_reference< decltype(quantityCallback.operator()(0))>::type
Calculates the value of an arbitrary vectorial quantity at any interior flux approximation point...
Definition: fvbasegradientcalculator.hh:146