AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
SecondOrderGradTestGradTrial.hpp
1#pragma once
2
3#include <type_traits>
4
5#include <amdis/GridFunctionOperator.hpp>
6#include <amdis/Output.hpp>
7#include <amdis/common/StaticSize.hpp>
8#include <amdis/common/ValueCategory.hpp>
9#include <amdis/typetree/FiniteElementType.hpp>
10
11namespace AMDiS
12{
18 namespace tag
19 {
21 }
22
23
26 {
27 public:
29
30 template <class CG, class RN, class CN, class Quad, class LocalFct, class Mat>
31 void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
32 Quad const& quad, LocalFct const& localFct, Mat& elementMatrix) const
33 {
34 using expr_value_type = typename LocalFct::Range;
35 static_assert(static_size_v<expr_value_type> == 1 ||
36 (static_num_rows_v<expr_value_type> == CG::dow &&
37 static_num_cols_v<expr_value_type> == CG::dow),
38 "Expression must be of scalar or matrix type.");
39 static_assert(RN::isLeaf && CN::isLeaf,
40 "Operator can be applied to Leaf-Nodes only.");
41
42 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
43 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
44 using Category = ValueCategory_t<typename LocalFct::Range>;
45
46 if (sameFE && sameNode)
47 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
48 else if (sameFE)
49 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
50 else
51 error_exit("Not implemented: currently only the implementation for equal fespaces available");
52 }
53
54 protected:
55
56 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
57 void getElementMatrixStandard(CG const& contextGeo, QR const& quad,
58 RN const& rowNode, CN const& colNode,
59 LocalFct const& localFct, Mat& elementMatrix) const
60 {
61 std::size_t size = rowNode.size();
62
63 using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
64 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
65 std::vector<WorldVector> gradients;
66
67 for (auto const& qp : quad) {
68 // Position of the current quadrature point in the reference element
69 auto&& local = contextGeo.coordinateInElement(qp.position());
70
71 // The transposed inverse Jacobian of the map from the reference element to the element
72 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
73
74 // The multiplicative factor in the integral transformation formula
75 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
76 const auto exprValue = localFct(local);
77
78 // The gradients of the shape functions on the reference element
79 auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
80
81 // Compute the shape function gradients on the real element
82 gradients.resize(shapeGradients.size());
83
84 for (std::size_t i = 0; i < gradients.size(); ++i)
85 jacobian.mv(shapeGradients[i][0], gradients[i]);
86
87 for (std::size_t i = 0; i < size; ++i) {
88 const auto local_i = rowNode.localIndex(i);
89 for (std::size_t j = 0; j < size; ++j) {
90 const auto local_j = colNode.localIndex(j);
91 elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
92 }
93 }
94 }
95 }
96
97 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
98 void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
99 RN const& node, CN const& /*colNode*/,
100 LocalFct const& localFct, Mat& elementMatrix, tag::scalar) const
101 {
102 std::size_t size = node.size();
103
104 using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
105 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
106 std::vector<WorldVector> gradients;
107
108 for (auto const& qp : quad) {
109 // Position of the current quadrature point in the reference element
110 auto&& local = contextGeo.coordinateInElement(qp.position());
111
112 // The transposed inverse Jacobian of the map from the reference element to the element
113 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
114
115 // The multiplicative factor in the integral transformation formula
116 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
117
118 // The gradients of the shape functions on the reference element
119 auto const& shapeGradients = node.localBasisJacobiansAt(local);
120
121 // Compute the shape function gradients on the real element
122 gradients.resize(shapeGradients.size());
123 for (std::size_t i = 0; i < gradients.size(); ++i)
124 jacobian.mv(shapeGradients[i][0], gradients[i]);
125
126 for (std::size_t i = 0; i < size; ++i) {
127 const auto local_i = node.localIndex(i);
128
129 elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
130
131 for (std::size_t j = i+1; j < size; ++j) {
132 const auto local_j = node.localIndex(j);
133 const auto value = factor * (gradients[i] * gradients[j]);
134
135 elementMatrix[local_i][local_j] += value;
136 elementMatrix[local_j][local_i] += value;
137 }
138 }
139 }
140 }
141
142 template <class CG, class QR, class RN, class CN, class LocalFct, class Mat>
143 void getElementMatrixOptimized(CG const& contextGeo, QR const& quad,
144 RN const& node, CN const& /*colNode*/,
145 LocalFct const& localFct, Mat& elementMatrix, tag::matrix) const
146 {
147 std::size_t size = node.size();
148
149 using RangeFieldType = typename RN::LocalBasis::Traits::RangeFieldType;
150 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
151 std::vector<WorldVector> gradients;
152
153 for (auto const& qp : quad) {
154 // Position of the current quadrature point in the reference element
155 auto&& local = contextGeo.coordinateInElement(qp.position());
156
157 // The transposed inverse Jacobian of the map from the reference element to the element
158 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
159
160 // The multiplicative factor in the integral transformation formula
161 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
162 const auto exprValue = localFct(local);
163
164 // The gradients of the shape functions on the reference element
165 auto const& shapeGradients = node.localBasisJacobiansAt(local);
166
167 // Compute the shape function gradients on the real element
168 gradients.resize(shapeGradients.size());
169 for (std::size_t i = 0; i < gradients.size(); ++i)
170 jacobian.mv(shapeGradients[i][0], gradients[i]);
171
172 for (std::size_t i = 0; i < size; ++i) {
173 const auto local_i = node.localIndex(i);
174 for (std::size_t j = 0; j < size; ++j) {
175 const auto local_j = node.localIndex(j);
176 elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
177 }
178 }
179 }
180 }
181
182 protected:
183
184 template <class S, class F, class T, int dow,
185 std::enable_if_t<Category::Scalar<S>,int> = 0>
186 T eval(S const& scalar, F factor,
187 Dune::FieldVector<T,dow> const& grad_test,
188 Dune::FieldVector<T,dow> const& grad_trial) const
189 {
190 return (scalar * factor) * (grad_test * grad_trial);
191 }
192
193 template <class M, class F, class T, int dow,
194 std::enable_if_t<Category::Matrix<M>,int> = 0>
195 T eval(M const& mat, F factor,
196 Dune::FieldVector<T,dow> const& grad_test,
197 Dune::FieldVector<T,dow> const& grad_trial) const
198 {
199 return factor * (grad_test * (mat * grad_trial));
200 }
201 };
202
203 template <class LC>
204 struct GridFunctionOperatorRegistry<tag::gradtest_gradtrial, LC>
205 {
206 static constexpr int degree = 2;
208 };
209
211 template <class Expr>
212 auto sot(Expr&& expr, int quadOrder = -1)
213 {
214 return makeOperator(tag::gradtest_gradtrial{}, FWD(expr), quadOrder);
215 }
216
219} // end namespace AMDiS
second-order operator , or
Definition: SecondOrderGradTestGradTrial.hpp:26
auto makeOperator(Tag const &tag, Expr &&expr, int gridFctDeg=-1)
Definition: GridFunctionOperator.hpp:235
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: SecondOrderGradTestGradTrial.hpp:20
Definition: Tags.hpp:11
Definition: Tags.hpp:9