AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
ZeroOrderTestTrial.hpp
1#pragma once
2
3#include <type_traits>
4
5#include <amdis/GridFunctionOperator.hpp>
6#include <amdis/common/StaticSize.hpp>
7#include <amdis/typetree/FiniteElementType.hpp>
8
9namespace AMDiS
10{
16 namespace tag
17 {
18 struct test_trial {};
19 }
20
21
24 {
25 public:
27
28 template <class CG, class RN, class CN, class Quad, class LocalFct, class Mat>
29 void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
30 Quad const& quad, LocalFct const& localFct, Mat& elementMatrix) const
31 {
32 static_assert(static_size_v<typename LocalFct::Range> == 1,
33 "Expression must be of scalar type." );
34
35 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
36 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
37
38 if (sameFE && sameNode)
39 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
40 else
41 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
42 }
43
44 template <class CG, class Node, class Quad, class LocalFct, class Vec>
45 void assemble(CG const& contextGeo, Node const& node,
46 Quad const& quad, LocalFct const& localFct, Vec& elementVector) const
47 {
48 static_assert(static_size_v<typename LocalFct::Range> == 1,
49 "Expression must be of scalar type." );
50 static_assert(Node::isLeaf,
51 "Operator can be applied to Leaf-Nodes only");
52
53 std::size_t size = node.size();
54
55 for (auto const& qp : quad) {
56 // Position of the current quadrature point in the reference element
57 auto&& local = contextGeo.coordinateInElement(qp.position());
58
59 // The multiplicative factor in the integral transformation formula
60 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
61
62 auto const& shapeValues = node.localBasisValuesAt(local);
63 for (std::size_t i = 0; i < size; ++i) {
64 const auto local_i = node.localIndex(i);
65 elementVector[local_i] += factor * shapeValues[i];
66 }
67 }
68 }
69
70
71 protected:
72 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
73 void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
74 RN const& rowNode, CN const& colNode,
75 LocalFct const& localFct, Mat& elementMatrix) const
76 {
77 static_assert(RN::isLeaf && CN::isLeaf,
78 "Operator can be applied to Leaf-Nodes only.");
79
80 std::size_t rowSize = rowNode.size();
81 std::size_t colSize = colNode.size();
82
83 for (auto const& qp : quad) {
84 // Position of the current quadrature point in the reference element
85 auto&& local = contextGeo.coordinateInElement(qp.position());
86
87 // The multiplicative factor in the integral transformation formula
88 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
89
90 auto const& rowShapeValues = rowNode.localBasisValuesAt(local);
91 auto const& colShapeValues = colNode.localBasisValuesAt(local);
92
93 for (std::size_t i = 0; i < rowSize; ++i) {
94 const auto local_i = rowNode.localIndex(i);
95 const auto value = factor * rowShapeValues[i];
96
97 for (std::size_t j = 0; j < colSize; ++j) {
98 const auto local_j = colNode.localIndex(j);
99 elementMatrix[local_i][local_j] += value * colShapeValues[j];
100 }
101 }
102 }
103 }
104
105
106 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
107 void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
108 RN const& node, CN const& /*colNode*/,
109 LocalFct const& localFct, Mat& elementMatrix) const
110 {
111 static_assert(RN::isLeaf && CN::isLeaf,
112 "Operator can be applied to Leaf-Nodes only.");
113
114 std::size_t size = node.size();
115
116 for (auto const& qp : quad) {
117 // Position of the current quadrature point in the reference element
118 auto&& local = contextGeo.coordinateInElement(qp.position());
119
120 // The multiplicative factor in the integral transformation formula
121 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
122
123 auto const& shapeValues = node.localBasisValuesAt(local);
124
125 for (std::size_t i = 0; i < size; ++i) {
126 const auto local_i = node.localIndex(i);
127
128 const auto value = factor * shapeValues[i];
129 elementMatrix[local_i][local_i] += value * shapeValues[i];
130
131 for (std::size_t j = i+1; j < size; ++j) {
132 const auto local_j = node.localIndex(j);
133
134 elementMatrix[local_i][local_j] += value * shapeValues[j];
135 elementMatrix[local_j][local_i] += value * shapeValues[j];
136 }
137 }
138 }
139 }
140 };
141
142 template <class LC>
143 struct GridFunctionOperatorRegistry<tag::test_trial, LC>
144 {
145 static constexpr int degree = 0;
146 using type = ZeroOrderTestTrial;
147 };
148
149
151 template <class Expr>
152 auto zot(Expr&& expr, int quadOrder = -1)
153 {
154 return makeOperator(tag::test_trial{}, FWD(expr), quadOrder);
155 }
156
159} // end namespace AMDiS
zero-order operator
Definition: ZeroOrderTestTrial.hpp:24
auto makeOperator(Tag const &tag, Expr &&expr, int gridFctDeg=-1)
Definition: GridFunctionOperator.hpp:235
auto zot(Expr &&expr, int quadOrder=-1)
Create a zero-order term.
Definition: ZeroOrderTestTrial.hpp:152
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: ZeroOrderTestTrial.hpp:18