5#include <dune/common/hybridutilities.hh>
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>
19 template <
class LocalFctA,
class LocalFctB,
class LocalFctC,
class LocalFctF,
20 bool conserving =
true>
21 class ConvectionDiffusionLocalOperator;
23 template <
class GridFctA,
class GridFctB,
class GridFctC,
class GridFctF,
24 bool conserving =
true>
30 GridFctC
const& gridFctC, GridFctF
const& gridFctF,
31 int quadOrder, bool_t<conserving> = {})
36 , quadOrder_(quadOrder)
39 template <
class Gr
idView>
40 void update(GridView
const&) { }
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>{}};
63 template <
class LocalFctA,
class LocalFctB,
class LocalFctC,
class LocalFctF,
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)
79 template <
class Element>
80 void bind(Element
const& element)
82 localFctA_.bind(element);
83 localFctB_.bind(element);
84 localFctC_.bind(element);
85 localFctF_.bind(element);
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
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." );
115 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
116 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
117 std::vector<WorldVector> gradients;
119 std::size_t numLocalFe = rowNode.size();
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)});
127 using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::Geometry::mydimension>;
128 auto const& quad = QuadratureRules::rule(contextGeo.type(), quadDeg);
130 for (std::size_t iq = 0; iq < quad.size(); ++iq) {
132 auto&& local = contextGeo.coordinateInElement(quad[iq].position());
135 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
138 const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
141 auto const& shapeValues = rowNode.localBasisValuesAt(local);
144 auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
147 gradients.resize(shapeGradients.size());
149 for (std::size_t i = 0; i < gradients.size(); ++i)
150 jacobian.mv(shapeGradients[i][0], gradients[i]);
152 const auto A = localFctA_(local);
153 WorldVector b = makeB<CG::dow>(localFctB_(local));
154 const auto c = localFctC_(local);
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;
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;
183 template <
class CG,
class Node,
class Vec>
184 void assemble(CG
const& contextGeo, Node
const& node, Vec& elementVector)
const
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." );
192 std::size_t numLocalFe = node.size();
194 auto const& quad = getQuadratureRule(contextGeo.geometry(), 0, coeffOrder(localFctF_), node);
196 for (std::size_t iq = 0; iq < quad.size(); ++iq) {
198 auto&& local = contextGeo.coordinateInElement(quad[iq].position());
201 auto const& shapeValues = node.localBasisValuesAt(local);
204 const auto factor = localFctF_(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
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;
214 template <
class LocalFct>
215 int coeffOrder(LocalFct
const& localFct)
const
217 if constexpr (Concepts::Polynomial<LocalFct>)
218 return order(localFct);
223 template <
int dow,
class T,
int N>
224 static FieldVector<T,dow> makeB(FieldVector<T,N>
const& b) {
return b; }
226 template <
int dow,
class T,
int N>
227 static FieldVector<T,dow> makeB(FieldVector<T,N>&& b) {
return std::move(b); }
229 template <
int dow,
class T>
230 static FieldVector<T,dow> makeB(FieldVector<T,1>
const& b) {
return FieldVector<T,dow>(T(b)); }
232 template <
int dow,
class T>
233 static FieldVector<T,dow> makeB(FieldVector<T,1>&& b) {
return FieldVector<T,dow>(T(b)); }
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); }
239 LocalFctA localFctA_;
240 LocalFctB localFctB_;
241 LocalFctC localFctC_;
242 LocalFctF localFctF_;
247 template <
class PreGr
idFctA,
class PreGr
idFctB,
class PreGr
idFctC,
class PreGr
idFctF,
class c>
250 static constexpr bool conserving = c::value;
251 PreGridFctA gridFctA;
252 PreGridFctB gridFctB;
253 PreGridFctC gridFctC;
254 PreGridFctF gridFctF;
258 template <
class Context,
class... T,
class GridView>
271 template <
bool conserving =
true,
272 class PreGridFctA,
class PreGridFctB,
class PreGridFctC,
class PreGridFctF>
274 PreGridFctC
const& gridFctC, PreGridFctF
const& gridFctF,
275 int quadOrder = -1, bool_t<conserving> = {})
277 using Pre = ConvectionDiffusionOperatorTerm<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, bool_t<conserving>>;
278 return Pre{gridFctA, gridFctB, gridFctC, gridFctF, quadOrder};
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