AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
ConvectionDiffusionOperator.hpp
1 #pragma once
2 
3 #include <type_traits>
4 
5 #include <dune/common/hybridutilities.hh>
6 
7 #include <amdis/LocalOperator.hpp>
8 #include <amdis/common/Order.hpp>
9 #include <amdis/common/StaticSize.hpp>
10 #include <amdis/common/TypeTraits.hpp>
11 #include <amdis/typetree/FiniteElementType.hpp>
12 
13 namespace AMDiS
14 {
19  template <class LocalFctA, class LocalFctB, class LocalFctC, class LocalFctF,
20  bool conserving = true>
22 
23  template <class GridFctA, class GridFctB, class GridFctC, class GridFctF,
24  bool conserving = true>
26  {
27  public:
29  ConvectionDiffusionOperator(GridFctA const& gridFctA, GridFctB const& gridFctB,
30  GridFctC const& gridFctC, GridFctF const& gridFctF,
31  int quadOrder, bool_t<conserving> = {})
32  : gridFctA_(gridFctA)
33  , gridFctB_(gridFctB)
34  , gridFctC_(gridFctC)
35  , gridFctF_(gridFctF)
36  , quadOrder_(quadOrder)
37  {}
38 
39  template <class GridView>
40  void update(GridView const&) { /* do nothing */ }
41 
42  friend auto localOperator(ConvectionDiffusionOperator const& op)
43  {
45  localFunction(op.gridFctA_),
46  localFunction(op.gridFctB_),
47  localFunction(op.gridFctC_),
48  localFunction(op.gridFctF_),
49  op.quadOrder_, bool_t<conserving>{}};
50  }
51 
52  private:
53  GridFctA gridFctA_;
54  GridFctB gridFctB_;
55  GridFctC gridFctC_;
56  GridFctF gridFctF_;
57  int quadOrder_;
58  };
59 
63  template <class LocalFctA, class LocalFctB, class LocalFctC, class LocalFctF,
64  bool conserving>
66  {
67  public:
69  LocalFctA const& localFctA, LocalFctB const& localFctB,
70  LocalFctC const& localFctC, LocalFctF const& localFctF,
71  int quadOrder, bool_t<conserving> = {})
72  : localFctA_(localFctA)
73  , localFctB_(localFctB)
74  , localFctC_(localFctC)
75  , localFctF_(localFctF)
76  , quadOrder_(quadOrder)
77  {}
78 
79  template <class Element>
80  void bind(Element const& element)
81  {
82  localFctA_.bind(element);
83  localFctB_.bind(element);
84  localFctC_.bind(element);
85  localFctF_.bind(element);
86  }
87 
89  void unbind()
90  {
91  localFctA_.unbind();
92  localFctB_.unbind();
93  localFctC_.unbind();
94  localFctF_.unbind();
95  }
96 
97  template <class CG, class RN, class CN, class Mat>
98  void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
99  Mat& elementMatrix) const
100  {
101  using A_range = typename LocalFctA::Range;
102  static_assert(static_size_v<A_range> == 1 || (static_num_rows_v<A_range> == CG::dow && static_num_cols_v<A_range> == CG::dow),
103  "Expression A(x) must be of scalar or matrix type.");
104  using b_range = typename LocalFctB::Range;
105  static_assert(static_size_v<b_range> == 1 || static_size_v<b_range> == CG::dow,
106  "Expression b(x) must be of scalar or vector type.");
107  using c_range = typename LocalFctC::Range;
108  static_assert(static_size_v<c_range> == 1,
109  "Expression c(x) must be of scalar type.");
110  static_assert(RN::isLeaf && CN::isLeaf,
111  "Operator can be applied to Leaf-Nodes only." );
112  static_assert(std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>,
113  "Galerkin operator requires equal finite elements for test and trial space." );
114 
115  using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
116  using WorldVector = FieldVector<RangeFieldType,CG::dow>;
117  std::vector<WorldVector> gradients;
118 
119  std::size_t numLocalFe = rowNode.size();
120 
121  auto contextGeometry = contextGeo.geometry();
122  int quadDeg = std::max({
123  getQuadratureDegree(contextGeometry,2,coeffOrder(localFctA_),rowNode,colNode),
124  getQuadratureDegree(contextGeometry,1,coeffOrder(localFctB_),rowNode,colNode),
125  getQuadratureDegree(contextGeometry,1,coeffOrder(localFctC_),rowNode,colNode)});
126 
127  using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::Geometry::mydimension>;
128  auto const& quad = QuadratureRules::rule(contextGeo.type(), quadDeg);
129 
130  for (std::size_t iq = 0; iq < quad.size(); ++iq) {
131  // Position of the current quadrature point in the reference element
132  auto&& local = contextGeo.coordinateInElement(quad[iq].position());
133 
134  // The transposed inverse Jacobian of the map from the reference element to the element
135  const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
136 
137  // The multiplicative factor in the integral transformation formula
138  const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
139 
140  // the values of the shape functions on the reference element at the quadrature point
141  auto const& shapeValues = rowNode.localBasisValuesAt(local);
142 
143  // The gradients of the shape functions on the reference element
144  auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
145 
146  // Compute the shape function gradients on the real element
147  gradients.resize(shapeGradients.size());
148 
149  for (std::size_t i = 0; i < gradients.size(); ++i)
150  jacobian.mv(shapeGradients[i][0], gradients[i]);
151 
152  const auto A = localFctA_(local);
153  WorldVector b = makeB<CG::dow>(localFctB_(local));
154  const auto c = localFctC_(local);
155 
156  if (conserving) {
157  WorldVector gradAi;
158  for (std::size_t i = 0; i < numLocalFe; ++i) {
159  const auto local_i = rowNode.localIndex(i);
160  gradAi = A * gradients[i];
161  auto gradBi = b * gradients[i];
162  for (std::size_t j = 0; j < numLocalFe; ++j) {
163  const auto local_j = colNode.localIndex(j);
164  elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
165  }
166  }
167  } else {
168  WorldVector grad_i;
169  for (std::size_t i = 0; i < numLocalFe; ++i) {
170  const auto local_i = rowNode.localIndex(i);
171  grad_i = A * gradients[i];
172  grad_i+= b * shapeValues[i];
173  for (std::size_t j = 0; j < numLocalFe; ++j) {
174  const auto local_j = colNode.localIndex(j);
175  elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
176  }
177  }
178  }
179  }
180  }
181 
182 
183  template <class CG, class Node, class Vec>
184  void assemble(CG const& contextGeo, Node const& node, Vec& elementVector) const
185  {
186  using f_range = typename LocalFctF::Range;
187  static_assert( static_size_v<f_range> == 1,
188  "Expression f(x) must be of scalar type." );
189  static_assert(Node::isLeaf,
190  "Operator can be applied to Leaf-Nodes only." );
191 
192  std::size_t numLocalFe = node.size();
193 
194  auto const& quad = getQuadratureRule(contextGeo.geometry(), 0, coeffOrder(localFctF_), node);
195 
196  for (std::size_t iq = 0; iq < quad.size(); ++iq) {
197  // Position of the current quadrature point in the reference element
198  auto&& local = contextGeo.coordinateInElement(quad[iq].position());
199 
200  // the values of the shape functions on the reference element at the quadrature point
201  auto const& shapeValues = node.localBasisValuesAt(local);
202 
203  // The multiplicative factor in the integral transformation formula
204  const auto factor = localFctF_(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
205 
206  for (std::size_t i = 0; i < numLocalFe; ++i) {
207  const auto local_i = node.localIndex(i);
208  elementVector[local_i] += shapeValues[i] * factor;
209  }
210  }
211  }
212 
213  private:
214  template <class LocalFct>
215  int coeffOrder(LocalFct const& localFct) const
216  {
217  if constexpr (Concepts::Polynomial<LocalFct>)
218  return order(localFct);
219  else
220  return 0;
221  }
222 
223  template <int dow, class T, int N>
224  static FieldVector<T,dow> makeB(FieldVector<T,N> const& b) { return b; }
225 
226  template <int dow, class T, int N>
227  static FieldVector<T,dow> makeB(FieldVector<T,N>&& b) { return std::move(b); }
228 
229  template <int dow, class T>
230  static FieldVector<T,dow> makeB(FieldVector<T,1> const& b) { return FieldVector<T,dow>(T(b)); }
231 
232  template <int dow, class T>
233  static FieldVector<T,dow> makeB(FieldVector<T,1>&& b) { return FieldVector<T,dow>(T(b)); }
234 
235  template <int dow, class T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
236  static FieldVector<T,dow> makeB(T b) { return FieldVector<T,dow>(b); }
237 
238  private:
239  LocalFctA localFctA_;
240  LocalFctB localFctB_;
241  LocalFctC localFctC_;
242  LocalFctF localFctF_;
243  int quadOrder_;
244  };
245 
246 #ifndef DOXYGEN
247  template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, class c>
249  {
250  static constexpr bool conserving = c::value;
251  PreGridFctA gridFctA;
252  PreGridFctB gridFctB;
253  PreGridFctC gridFctC;
254  PreGridFctF gridFctF;
255  int quadOrder;
256  };
257 
258  template <class Context, class... T, class GridView>
259  auto makeOperator(ConvectionDiffusionOperatorTerm<T...> const& pre, GridView const& gridView)
260  {
262  makeGridFunction(pre.gridFctA, gridView),
263  makeGridFunction(pre.gridFctB, gridView),
264  makeGridFunction(pre.gridFctC, gridView),
265  makeGridFunction(pre.gridFctF, gridView),
266  pre.quadOrder};
267  }
268 #endif
269 
271  template <bool conserving = true,
272  class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF>
273  auto convectionDiffusion(PreGridFctA const& gridFctA, PreGridFctB const& gridFctB,
274  PreGridFctC const& gridFctC, PreGridFctF const& gridFctF,
275  int quadOrder = -1, bool_t<conserving> = {})
276  {
278  return Pre{gridFctA, gridFctB, gridFctC, gridFctF, quadOrder};
279  }
280 
283 } // end namespace AMDiS
ConvectionDiffusionOperator(GridFctA const &gridFctA, GridFctB const &gridFctB, GridFctC const &gridFctC, GridFctF const &gridFctF, int quadOrder, bool_t< conserving >={})
Constructor. Stores a copy of all passed gridfunctions.
Definition: ConvectionDiffusionOperator.hpp:29
decltype(auto) makeGridFunction(PreGridFct const &preGridFct, GridView const &gridView)
Generator for Gridfunctions from Expressions (PreGridfunctions)
Definition: GridFunction.hpp:168
Definition: AdaptBase.hpp:6
Definition: ConvectionDiffusionOperator.hpp:21
void unbind()
Unbinds operator from element.
Definition: ConvectionDiffusionOperator.hpp:89
auto makeOperator(Tag const &tag, Expr &&expr, int gridFctDeg=-1)
Definition: GridFunctionOperator.hpp:235
Definition: ConvectionDiffusionOperator.hpp:248
auto convectionDiffusion(PreGridFctA const &gridFctA, PreGridFctB const &gridFctB, PreGridFctC const &gridFctC, PreGridFctF const &gridFctF, int quadOrder=-1, bool_t< conserving >={})
Generator function for constructing a Convection-Diffusion Operator.
Definition: ConvectionDiffusionOperator.hpp:273
Definition: ConvectionDiffusionOperator.hpp:25