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>
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
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.");
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>;
46 if (sameFE && sameNode)
47 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
49 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
51 error_exit(
"Not implemented: currently only the implementation for equal fespaces available");
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
61 std::size_t size = rowNode.size();
63 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
64 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
65 std::vector<WorldVector> gradients;
67 for (
auto const& qp : quad) {
69 auto&& local = contextGeo.coordinateInElement(qp.position());
72 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
75 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
76 const auto exprValue = localFct(local);
79 auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
82 gradients.resize(shapeGradients.size());
84 for (std::size_t i = 0; i < gradients.size(); ++i)
85 jacobian.mv(shapeGradients[i][0], gradients[i]);
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]);
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& ,
100 LocalFct
const& localFct, Mat& elementMatrix,
tag::scalar)
const
102 std::size_t size = node.size();
104 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
105 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
106 std::vector<WorldVector> gradients;
108 for (
auto const& qp : quad) {
110 auto&& local = contextGeo.coordinateInElement(qp.position());
113 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
116 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
119 auto const& shapeGradients = node.localBasisJacobiansAt(local);
122 gradients.resize(shapeGradients.size());
123 for (std::size_t i = 0; i < gradients.size(); ++i)
124 jacobian.mv(shapeGradients[i][0], gradients[i]);
126 for (std::size_t i = 0; i < size; ++i) {
127 const auto local_i = node.localIndex(i);
129 elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
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]);
135 elementMatrix[local_i][local_j] += value;
136 elementMatrix[local_j][local_i] += value;
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& ,
145 LocalFct
const& localFct, Mat& elementMatrix,
tag::matrix)
const
147 std::size_t size = node.size();
149 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
150 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
151 std::vector<WorldVector> gradients;
153 for (
auto const& qp : quad) {
155 auto&& local = contextGeo.coordinateInElement(qp.position());
158 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
161 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
162 const auto exprValue = localFct(local);
165 auto const& shapeGradients = node.localBasisJacobiansAt(local);
168 gradients.resize(shapeGradients.size());
169 for (std::size_t i = 0; i < gradients.size(); ++i)
170 jacobian.mv(shapeGradients[i][0], gradients[i]);
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]);
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
190 return (scalar * factor) * (grad_test * grad_trial);
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
199 return factor * (grad_test * (mat * grad_trial));
206 static constexpr int degree = 2;
211 template <
class Expr>
212 auto sot(Expr&& expr,
int quadOrder = -1)
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