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
13namespace AMDiS
14{
19 template <class LocalFctA, class LocalFctB, class LocalFctC, class LocalFctF,
20 bool conserving = true>
21 class ConvectionDiffusionLocalOperator;
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 {
44 return ConvectionDiffusionLocalOperator{
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 quadOrder_;
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 {
277 using Pre = ConvectionDiffusionOperatorTerm<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, bool_t<conserving>>;
278 return Pre{gridFctA, gridFctB, gridFctC, gridFctF, quadOrder};
279 }
280
283} // end namespace AMDiS
Definition: ConvectionDiffusionOperator.hpp:66
void unbind()
Unbinds operator from element.
Definition: ConvectionDiffusionOperator.hpp:89
Definition: ConvectionDiffusionOperator.hpp:26
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
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
auto makeOperator(Tag const &tag, Expr &&expr, int gridFctDeg=-1)
Definition: GridFunctionOperator.hpp:235
Definition: ConvectionDiffusionOperator.hpp:249