AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
ZeroOrderTestvecTrialvec.hpp
1#pragma once
2
3#include <type_traits>
4
5#include <amdis/GridFunctionOperator.hpp>
6#include <amdis/common/StaticSize.hpp>
7#include <amdis/common/ValueCategory.hpp>
8#include <amdis/typetree/FiniteElementType.hpp>
9
10namespace AMDiS
11{
17 namespace tag
18 {
20 }
21
22
25 {
26 public:
28
29 template <class CG, class RN, class CN, class Quad, class LocalFct, class Mat>
30 void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
31 Quad const& quad, LocalFct const& localFct, Mat& elementMatrix) const
32 {
33 using expr_value_type = typename LocalFct::Range;
34 static_assert(static_size_v<expr_value_type> == 1 ||
35 (static_num_rows_v<expr_value_type> == CG::dow &&
36 static_num_cols_v<expr_value_type> == CG::dow),
37 "Expression must be of scalar or matrix type." );
38 static_assert(RN::isPower && CN::isPower,
39 "Operator can be applied to Power-Nodes only.");
40 assert(rowNode.degree() == colNode.degree());
41
42 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
43 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
44
45 using Category = ValueCategory_t<expr_value_type>;
46
47 if (sameFE && sameNode && std::is_same_v<Category,tag::scalar>)
48 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
49 else
50 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
51 }
52
53 protected: // specialization for different value-categories
54
55 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
56 void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
57 RN const& rowNode, CN const& colNode,
58 LocalFct const& localFct, Mat& elementMatrix, tag::scalar) const
59 {
60 std::size_t rowSize = rowNode.child(0).size();
61 std::size_t colSize = colNode.child(0).size();
62
63 for (auto const& qp : quad) {
64 // Position of the current quadrature point in the reference element
65 auto&& local = contextGeo.coordinateInElement(qp.position());
66
67 // The multiplicative factor in the integral transformation formula
68 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
69
70 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
71 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
72
73 for (std::size_t i = 0; i < rowSize; ++i) {
74 const auto value = factor * rowShapeValues[i];
75
76 for (std::size_t j = 0; j < colSize; ++j) {
77 const auto value0 = value * colShapeValues[j];
78
79 for (std::size_t k = 0; k < rowNode.degree(); ++k) {
80 const auto local_ki = rowNode.child(k).localIndex(i);
81 const auto local_kj = colNode.child(k).localIndex(j);
82
83 elementMatrix[local_ki][local_kj] += value0;
84 }
85 }
86 }
87 }
88 }
89
90
91 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
92 void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
93 RN const& node, CN const& /*colNode*/,
94 LocalFct const& localFct, Mat& elementMatrix, tag::scalar) const
95 {
96 std::size_t size = node.child(0).size();
97
98 for (auto const& qp : quad) {
99 // Position of the current quadrature point in the reference element
100 auto&& local = contextGeo.coordinateInElement(qp.position());
101
102 // The multiplicative factor in the integral transformation formula
103 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
104
105 auto const& shapeValues = node.child(0).localBasisValuesAt(local);
106
107 for (std::size_t i = 0; i < size; ++i) {
108 const auto value = factor * shapeValues[i];
109
110 for (std::size_t k = 0; k < node.degree(); ++k) {
111 const auto local_ki = node.child(k).localIndex(i);
112 elementMatrix[local_ki][local_ki] += value * shapeValues[i];
113 }
114
115 for (std::size_t j = i+1; j < size; ++j) {
116 const auto value0 = value * shapeValues[j];
117
118 for (std::size_t k = 0; k < node.degree(); ++k) {
119 const auto local_ki = node.child(k).localIndex(i);
120 const auto local_kj = node.child(k).localIndex(j);
121
122 elementMatrix[local_ki][local_kj] += value0;
123 elementMatrix[local_kj][local_ki] += value0;
124 }
125 }
126 }
127 }
128 }
129
130 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
131 void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
132 RN const& rowNode, CN const& colNode,
133 LocalFct const& localFct, Mat& elementMatrix, tag::matrix) const
134 {
135 using expr_value_type [[maybe_unused]] = typename LocalFct::Range;
136 assert(static_num_rows_v<expr_value_type> == rowNode.degree() &&
137 static_num_cols_v<expr_value_type> == rowNode.degree());
138
139 std::size_t rowSize = rowNode.child(0).size();
140 std::size_t colSize = colNode.child(0).size();
141
142 for (auto const& qp : quad) {
143 // Position of the current quadrature point in the reference element
144 auto&& local = contextGeo.coordinateInElement(qp.position());
145
146 // The multiplicative factor in the integral transformation formula
147 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
148 const auto exprValue = localFct(local);
149
150 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
151 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
152
153 for (std::size_t i = 0; i < rowSize; ++i) {
154 const auto value0 = factor * rowShapeValues[i];
155
156 for (std::size_t j = 0; j < colSize; ++j) {
157 const auto value = value0 * colShapeValues[j];
158 const auto mat = exprValue * value;
159
160 for (std::size_t k0 = 0; k0 < rowNode.degree(); ++k0) {
161 const auto local_ki = rowNode.child(k0).localIndex(i);
162 for (std::size_t k1 = 0; k1 < rowNode.degree(); ++k1) {
163 const auto local_kj = colNode.child(k1).localIndex(j);
164
165 elementMatrix[local_ki][local_kj] += mat[k0][k1];
166 }
167 }
168 }
169 }
170 }
171 }
172
173 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
174 void getElementMatrixOptimized(CG const&, QR const&, RN const&, CN const&,
175 LocalFct const&, Mat&, tag::matrix) const
176 {
177 DUNE_THROW(Dune::NotImplemented,
178 "Optimized LocalOperator with matrix coefficients not implemented.");
179 }
180 };
181
182 template <class LC>
183 struct GridFunctionOperatorRegistry<tag::testvec_trialvec, LC>
184 {
185 static constexpr int degree = 0;
187 };
188
191} // end namespace AMDiS
zero-order operator , or
Definition: ZeroOrderTestvecTrialvec.hpp:25
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: Tags.hpp:11
Definition: Tags.hpp:9
Definition: ZeroOrderTestvecTrialvec.hpp:19