5 #include <dune/common/hybridutilities.hh> 7 #include <amdis/LocalOperator.hpp> 8 #include <amdis/Output.hpp> 9 #include <amdis/common/TypeTraits.hpp> 10 #include <amdis/common/StaticSize.hpp> 19 template <
class LC,
class Gr
idFctA,
class Gr
idFctB,
class Gr
idFctC,
class Gr
idFctF,
bool conserving = true>
24 :
public LocalOperator<ConvectionDiffusionOperator<LC, GridFctA, GridFctB, GridFctC, GridFctF, conserving>, LC>
28 using A_range =
typename GridFctA::Range;
29 static_assert( static_size_v<A_range> == 1 || (static_num_rows_v<A_range> == dow && static_num_cols_v<A_range> == dow),
30 "Expression A(x) must be of scalar or matrix type." );
31 using b_range =
typename GridFctB::Range;
32 static_assert( static_size_v<b_range> == 1 || static_size_v<b_range> == dow,
33 "Expression b(x) must be of scalar or vector type." );
34 using c_range =
typename GridFctC::Range;
35 static_assert( static_size_v<c_range> == 1,
36 "Expression c(x) must be of scalar type." );
37 using f_range =
typename GridFctF::Range;
38 static_assert( static_size_v<f_range> == 1,
39 "Expression f(x) must be of scalar type." );
43 GridFctC
const& gridFctC, GridFctF
const& gridFctF)
50 template <
class CG,
class RN,
class CN,
class Mat>
51 void getElementMatrix(CG
const& contextGeo, RN
const& rowNode, CN
const& colNode, Mat& elementMatrix)
53 static_assert(RN::isLeaf && CN::isLeaf,
54 "Operator can be applied to Leaf-Nodes only." );
56 static_assert(std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>,
57 "Galerkin operator requires equal finite elements for test and trial space." );
59 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
60 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
61 std::vector<WorldVector> gradients;
63 auto localFctA = localFunction(gridFctA_); localFctA.bind(contextGeo.element());
64 auto localFctB = localFunction(gridFctB_); localFctB.bind(contextGeo.element());
65 auto localFctC = localFunction(gridFctC_); localFctC.bind(contextGeo.element());
67 std::size_t numLocalFe = rowNode.size();
69 int quadDeg = std::max({this->
getDegree(2,coeffOrder(localFctA),rowNode,colNode),
70 this->
getDegree(1,coeffOrder(localFctB),rowNode,colNode),
71 this->
getDegree(0,coeffOrder(localFctC),rowNode,colNode)});
73 using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::LocalContext::mydimension>;
74 auto const& quad = QuadratureRules::rule(contextGeo.type(), quadDeg);
76 for (std::size_t iq = 0; iq < quad.size(); ++iq) {
78 auto&& local = contextGeo.local(quad[iq].position());
81 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
84 const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
87 auto const& shapeValues = rowNode.localBasisValuesAt(local);
90 auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
93 gradients.resize(shapeGradients.size());
95 for (std::size_t i = 0; i < gradients.size(); ++i)
96 jacobian.mv(shapeGradients[i][0], gradients[i]);
98 const auto A = localFctA(local);
99 WorldVector b = makeB(localFctB(local));
100 const auto c = localFctC(local);
104 for (std::size_t i = 0; i < numLocalFe; ++i) {
105 const auto local_i = rowNode.localIndex(i);
106 gradAi = A * gradients[i];
107 auto gradBi = b * gradients[i];
108 for (std::size_t j = 0; j < numLocalFe; ++j) {
109 const auto local_j = colNode.localIndex(j);
110 elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
115 for (std::size_t i = 0; i < numLocalFe; ++i) {
116 const auto local_i = rowNode.localIndex(i);
117 grad_i = A * gradients[i];
118 grad_i+= b * shapeValues[i];
119 for (std::size_t j = 0; j < numLocalFe; ++j) {
120 const auto local_j = colNode.localIndex(j);
121 elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
133 template <
class CG,
class Node,
class Vec>
134 void getElementVector(CG
const& contextGeo, Node
const& node, Vec& elementVector)
136 static_assert(Node::isLeaf,
137 "Operator can be applied to Leaf-Nodes only." );
139 auto localFctF = localFunction(gridFctF_); localFctF.bind(contextGeo.element());
141 std::size_t numLocalFe = node.size();
143 int quad_order = this->
getDegree(0,coeffOrder(localFctF),node);
145 using QuadratureRules = Dune::QuadratureRules<typename CG::Geometry::ctype, CG::LocalContext::dimension>;
146 auto const& quad = QuadratureRules::rule(contextGeo.type(), quad_order);
148 for (std::size_t iq = 0; iq < quad.size(); ++iq) {
150 auto&& local = contextGeo.local(quad[iq].position());
153 auto const& shapeValues = node.localBasisValuesAt(local);
156 const auto factor = localFctF(local) * contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
158 for (std::size_t i = 0; i < numLocalFe; ++i) {
159 const auto local_i = node.localIndex(i);
160 elementVector[local_i] += shapeValues[i] * factor;
169 template <
class LocalFct>
170 int coeffOrder(LocalFct
const& localFct)
172 if constexpr (Concepts::Polynomial<LocalFct>)
173 return order(localFct);
178 template <
class T,
int N>
179 static FieldVector<T,dow> makeB(FieldVector<T,N>
const& b) {
return b; }
181 template <
class T,
int N>
182 static FieldVector<T,dow> makeB(FieldVector<T,N>&& b) {
return std::move(b); }
185 static FieldVector<T,dow> makeB(FieldVector<T,1>
const& b) {
return FieldVector<T,dow>(T(b)); }
188 static FieldVector<T,dow> makeB(FieldVector<T,1>&& b) {
return FieldVector<T,dow>(T(b)); }
190 template <
class T, std::enable_if_t<std::is_arithmetic_v<T>,
int> = 0>
191 static FieldVector<T,dow> makeB(T b) {
return FieldVector<T,dow>(b); }
201 template <
class PreGr
idFctA,
class PreGr
idFctB,
class PreGr
idFctC,
class PreGr
idFctF,
class c>
204 static constexpr
bool conserving = c::value;
205 PreGridFctA gridFctA;
206 PreGridFctB gridFctB;
207 PreGridFctC gridFctC;
208 PreGridFctF gridFctF;
213 template <
class PreGr
idFctA,
class PreGr
idFctB,
class PreGr
idFctC,
class PreGr
idFctF,
bool conserving = true>
214 auto convectionDiffusion(PreGridFctA
const& gridFctA, PreGridFctB
const& gridFctB,
215 PreGridFctC
const& gridFctC, PreGridFctF
const& gridFctF,
219 return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
224 template <
class Context,
class... T,
class GridView>
234 decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), Pre::conserving>;
236 GridFctOp localOperator{std::move(gridFctA), std::move(gridFctB), std::move(gridFctC), std::move(gridFctF)};
237 return localOperator;
Definition: ConvectionDiffusionOperator.hpp:202
decltype(auto) makeGridFunction(PreGridFct const &preGridFct, GridView const &gridView)
Generator for Gridfunctions from Expressions (PreGridfunctions)
Definition: GridFunction.hpp:154
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
int getDegree(int derivOrder, int coeffDegree, RN const &rowNode, CN const &colNode) const
Return the quadrature degree for a matrix operator.
Definition: LocalOperator.hpp:130
auto order(F const &f) -> decltype(&F::operator(), f.order())
polynomial order of functions
Definition: Order.hpp:11
auto makeLocalOperator(PreGridFunctionOperator< Args... > op, GridView const &gridView)
Definition: GridFunctionOperator.hpp:258
The main implementation of an operator to be used in a Assembler.
Definition: LocalOperator.hpp:28
Wrapper class for element and geometry.
Definition: ContextGeometry.hpp:43
std::integral_constant< bool, B > bool_t
A wrapper for bool types.
Definition: Logical.hpp:12
Definition: ConvectionDiffusionOperator.hpp:23