All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
p1fegradientcalculator.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_P1FE_GRADIENT_CALCULATOR_HH
29 #define EWOMS_P1FE_GRADIENT_CALCULATOR_HH
30 
31 #include "vcfvproperties.hh"
32 
34 
35 #include <dune/common/fvector.hh>
36 #include <dune/common/version.hh>
37 
38 #include <dune/geometry/type.hh>
39 
40 #if HAVE_DUNE_LOCALFUNCTIONS
41 #include <dune/localfunctions/lagrange/pqkfactory.hh>
42 #endif // HAVE_DUNE_LOCALFUNCTIONS
43 
44 #include <dune/common/fvector.hh>
45 
46 #include <vector>
47 
48 #ifdef HAVE_DUNE_LOCALFUNCTIONS
49 #define EWOMS_NO_LOCALFUNCTIONS_UNUSED
50 #else
51 #define EWOMS_NO_LOCALFUNCTIONS_UNUSED OPM_UNUSED
52 #endif
53 
54 namespace Ewoms {
64 template<class TypeTag>
66 {
68  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
69  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
70  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
71  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
72 
73  enum { dim = GridView::dimension };
74 
75  // set the maximum number of degrees of freedom and the maximum
76  // number of flux approximation points per elements. For this, we
77  // assume cubes as the type of element with the most vertices...
78  enum { maxDof = (2 << dim) };
79  enum { maxFap = maxDof };
80 
81  typedef typename GridView::ctype CoordScalar;
82  typedef Dune::FieldVector<Scalar, dim> DimVector;
83 
84 #if HAVE_DUNE_LOCALFUNCTIONS
85  typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
86  typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
87  typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
88  typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
89 #endif // HAVE_DUNE_LOCALFUNCTIONS
90 
91 public:
98  template <bool prepareValues = true, bool prepareGradients = true>
99  void prepare(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
100  unsigned timeIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED)
101  {
102  if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
103 #if !HAVE_DUNE_LOCALFUNCTIONS
104  // The dune-localfunctions module is required for P1 finite element gradients
105  OPM_THROW(std::logic_error, "The dune-localfunctions module is required in oder to use"
106  " finite element gradients");
107 #else
108  const auto& stencil = elemCtx.stencil(timeIdx);
109 
110  const LocalFiniteElement& localFE = feCache_.get(elemCtx.element().type());
111  localFiniteElement_ = &localFE;
112 
113  // loop over all face centeres
114  for (unsigned faceIdx = 0; faceIdx < stencil.numInteriorFaces(); ++faceIdx) {
115  const auto& localFacePos = stencil.interiorFace(faceIdx).localPos();
116 
117  // Evaluate the P1 shape functions and their gradients at all
118  // flux approximation points.
119  if (prepareValues)
120  localFE.localBasis().evaluateFunction(localFacePos, p1Value_[faceIdx]);
121 
122  if (prepareGradients) {
123  // first, get the shape function's gradient in local coordinates
124  std::vector<ShapeJacobian> localGradient;
125  localFE.localBasis().evaluateJacobian(localFacePos, localGradient);
126 
127  // convert to a gradient in global space by
128  // multiplying with the inverse transposed jacobian of
129  // the position
130  const auto& geom = elemCtx.element().geometry();
131  const auto& jacInvT =
132  geom.jacobianInverseTransposed(localFacePos);
133 
134  size_t numVertices = elemCtx.numDof(timeIdx);
135  for (unsigned vertIdx = 0; vertIdx < numVertices; vertIdx++) {
136  jacInvT.mv(/*xVector=*/localGradient[vertIdx][0],
137  /*destVector=*/p1Gradient_[faceIdx][vertIdx]);
138  }
139  }
140  }
141 #endif
142  }
143  else
144  ParentType::template prepare<prepareValues, prepareGradients>(elemCtx, timeIdx);
145  }
146 
158  template <class QuantityCallback>
159  auto calculateScalarValue(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
160  unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
161  const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
162  -> typename std::remove_reference<typename QuantityCallback::ResultType>::type
163  {
164  if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
165 #if !HAVE_DUNE_LOCALFUNCTIONS
166  // The dune-localfunctions module is required for P1 finite element gradients
167  OPM_THROW(std::logic_error, "The dune-localfunctions module is required in oder to use"
168  " finite element gradients");
169 #else
170  typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
171  typedef typename std::remove_const<QuantityConstType>::type QuantityType;
172  typedef Opm::MathToolbox<QuantityType> Toolbox;
173 
174  // If the user does not want to use two-point gradients, we
175  // use P1 finite element gradients..
176  QuantityType value(0.0);
177  for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
178  if (std::is_same<QuantityType, Scalar>::value ||
179  elemCtx.focusDofIndex() == vertIdx)
180  value += quantityCallback(vertIdx)*p1Value_[fapIdx][vertIdx];
181  else
182  value += Toolbox::value(quantityCallback(vertIdx))*p1Value_[fapIdx][vertIdx];
183  }
184 
185  return value;
186 #endif
187  }
188  else
189  return ParentType::calculateScalarValue(elemCtx, fapIdx, quantityCallback);
190  }
191 
203  template <class QuantityCallback>
204  auto calculateVectorValue(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
205  unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
206  const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
207  -> typename std::remove_reference<typename QuantityCallback::ResultType>::type
208  {
209  if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
210 #if !HAVE_DUNE_LOCALFUNCTIONS
211  // The dune-localfunctions module is required for P1 finite element gradients
212  OPM_THROW(std::logic_error, "The dune-localfunctions module is required in oder to use"
213  " finite element gradients");
214 #else
215  typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
216  typedef typename std::remove_const<QuantityConstType>::type QuantityType;
217 
218  typedef decltype(std::declval<QuantityType>()[0]) RawFieldType;
219  typedef typename std::remove_const<typename std::remove_reference<RawFieldType>::type>::type FieldType;
220 
221  typedef Opm::MathToolbox<FieldType> Toolbox;
222 
223  // If the user does not want to use two-point gradients, we
224  // use P1 finite element gradients..
225  QuantityType value(0.0);
226  for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
227  if (std::is_same<QuantityType, Scalar>::value ||
228  elemCtx.focusDofIndex() == vertIdx)
229  {
230  const auto& tmp = quantityCallback(vertIdx);
231  for (unsigned k = 0; k < tmp.size(); ++k)
232  value[k] += tmp[k]*p1Value_[fapIdx][vertIdx];
233  }
234  else {
235  const auto& tmp = quantityCallback(vertIdx);
236  for (unsigned k = 0; k < tmp.size(); ++k)
237  value[k] += Toolbox::value(tmp[k])*p1Value_[fapIdx][vertIdx];
238  }
239  }
240 
241  return value;
242 #endif
243  }
244  else
245  return ParentType::calculateVectorValue(elemCtx, fapIdx, quantityCallback);
246  }
247 
259  template <class QuantityCallback, class EvalDimVector>
260  void calculateGradient(EvalDimVector& quantityGrad EWOMS_NO_LOCALFUNCTIONS_UNUSED,
261  const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
262  unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
263  const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
264  {
265  if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
266 #if !HAVE_DUNE_LOCALFUNCTIONS
267  // The dune-localfunctions module is required for P1 finite element gradients
268  OPM_THROW(std::logic_error, "The dune-localfunctions module is required in oder to use"
269  " finite element gradients");
270 #else
271  typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
272  typedef typename std::remove_const<QuantityConstType>::type QuantityType;
273 
274  // If the user does not want two-point gradients, we use P1 finite element
275  // gradients...
276  quantityGrad = 0.0;
277  for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
278  if (std::is_same<QuantityType, Scalar>::value ||
279  elemCtx.focusDofIndex() == vertIdx)
280  {
281  const auto& dofVal = quantityCallback(vertIdx);
282  const auto& tmp = p1Gradient_[fapIdx][vertIdx];
283  for (int dimIdx = 0; dimIdx < dim; ++ dimIdx)
284  quantityGrad[dimIdx] += dofVal*tmp[dimIdx];
285  }
286  else {
287  const auto& dofVal = quantityCallback(vertIdx);
288  const auto& tmp = p1Gradient_[fapIdx][vertIdx];
289  for (int dimIdx = 0; dimIdx < dim; ++ dimIdx)
290  quantityGrad[dimIdx] += Opm::scalarValue(dofVal)*tmp[dimIdx];
291  }
292  }
293 #endif
294  }
295  else
296  ParentType::calculateGradient(quantityGrad, elemCtx, fapIdx, quantityCallback);
297  }
298 
313  template <class QuantityCallback>
314  auto calculateBoundaryValue(const ElementContext& elemCtx,
315  unsigned fapIdx,
316  const QuantityCallback& quantityCallback)
317  -> decltype(ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback))
318  { return ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback); }
319 
334  template <class QuantityCallback, class EvalDimVector>
335  void calculateBoundaryGradient(EvalDimVector& quantityGrad,
336  const ElementContext& elemCtx,
337  unsigned fapIdx,
338  const QuantityCallback& quantityCallback) const
339  { ParentType::calculateBoundaryGradient(quantityGrad, elemCtx, fapIdx, quantityCallback); }
340 
341 #if HAVE_DUNE_LOCALFUNCTIONS
342  static LocalFiniteElementCache& localFiniteElementCache()
343  { return feCache_; }
344 #endif
345 
346 private:
347 #if HAVE_DUNE_LOCALFUNCTIONS
348  static LocalFiniteElementCache feCache_;
349 
350  const LocalFiniteElement* localFiniteElement_;
351  std::vector<Dune::FieldVector<Scalar, 1>> p1Value_[maxFap];
352  DimVector p1Gradient_[maxFap][maxDof];
353 #endif // HAVE_DUNE_LOCALFUNCTIONS
354 };
355 
356 #if HAVE_DUNE_LOCALFUNCTIONS
357 template<class TypeTag>
358 typename P1FeGradientCalculator<TypeTag>::LocalFiniteElementCache
359 P1FeGradientCalculator<TypeTag>::feCache_;
360 #endif
361 } // namespace Ewoms
362 
363 #undef EWOMS_NO_LOCALFUNCTIONS_UNUSED
364 
365 #endif
auto calculateVectorValue(const ElementContext &elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED, unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED, const QuantityCallback &quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const -> typename std::remove_reference< typename QuantityCallback::ResultType >::type
Calculates the value of an arbitrary quantity at any interior flux approximation point.
Definition: p1fegradientcalculator.hh:204
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
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
auto calculateBoundaryValue(const ElementContext &elemCtx, unsigned 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: p1fegradientcalculator.hh:314
void prepare(const ElementContext &elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED, unsigned timeIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED)
Precomputes the common values to calculate gradients and values of quantities at any flux approximati...
Definition: p1fegradientcalculator.hh:99
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 using the two-poin...
Definition: fvbasegradientcalculator.hh:48
This class calculates gradients of arbitrary quantities at flux integration points using first order ...
Definition: p1fegradientcalculator.hh:65
void calculateGradient(EvalDimVector &quantityGrad EWOMS_NO_LOCALFUNCTIONS_UNUSED, const ElementContext &elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED, unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED, const QuantityCallback &quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
Calculates the gradient of an arbitrary quantity at any flux approximation point. ...
Definition: p1fegradientcalculator.hh:260
Declares the basic properties used by the common infrastructure of the vertex-centered finite volume ...
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
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
void calculateBoundaryGradient(EvalDimVector &quantityGrad, const ElementContext &elemCtx, unsigned fapIdx, const QuantityCallback &quantityCallback) const
Calculates the gradient of an arbitrary quantity at any flux approximation point on the boundary...
Definition: p1fegradientcalculator.hh:335
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
auto calculateScalarValue(const ElementContext &elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED, unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED, const QuantityCallback &quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const -> typename std::remove_reference< typename QuantityCallback::ResultType >::type
Calculates the value of an arbitrary quantity at any interior flux approximation point.
Definition: p1fegradientcalculator.hh:159