5#include <amdis/GridFunctionOperator.hpp>
6#include <amdis/common/StaticSize.hpp>
7#include <amdis/common/ValueCategory.hpp>
8#include <amdis/typetree/FiniteElementType.hpp>
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
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());
42 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
43 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
45 using Category = ValueCategory_t<expr_value_type>;
47 if (sameFE && sameNode && std::is_same_v<Category,tag::scalar>)
48 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
50 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix, Category{});
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
60 std::size_t rowSize = rowNode.child(0).size();
61 std::size_t colSize = colNode.child(0).size();
63 for (
auto const& qp : quad) {
65 auto&& local = contextGeo.coordinateInElement(qp.position());
68 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
70 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
71 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
73 for (std::size_t i = 0; i < rowSize; ++i) {
74 const auto value = factor * rowShapeValues[i];
76 for (std::size_t j = 0; j < colSize; ++j) {
77 const auto value0 = value * colShapeValues[j];
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);
83 elementMatrix[local_ki][local_kj] += value0;
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& ,
94 LocalFct
const& localFct, Mat& elementMatrix,
tag::scalar)
const
96 std::size_t size = node.child(0).size();
98 for (
auto const& qp : quad) {
100 auto&& local = contextGeo.coordinateInElement(qp.position());
103 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
105 auto const& shapeValues = node.child(0).localBasisValuesAt(local);
107 for (std::size_t i = 0; i < size; ++i) {
108 const auto value = factor * shapeValues[i];
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];
115 for (std::size_t j = i+1; j < size; ++j) {
116 const auto value0 = value * shapeValues[j];
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);
122 elementMatrix[local_ki][local_kj] += value0;
123 elementMatrix[local_kj][local_ki] += value0;
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
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());
139 std::size_t rowSize = rowNode.child(0).size();
140 std::size_t colSize = colNode.child(0).size();
142 for (
auto const& qp : quad) {
144 auto&& local = contextGeo.coordinateInElement(qp.position());
147 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
148 const auto exprValue = localFct(local);
150 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
151 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
153 for (std::size_t i = 0; i < rowSize; ++i) {
154 const auto value0 = factor * rowShapeValues[i];
156 for (std::size_t j = 0; j < colSize; ++j) {
157 const auto value = value0 * colShapeValues[j];
158 const auto mat = exprValue * value;
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);
165 elementMatrix[local_ki][local_kj] += mat[k0][k1];
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&,
177 DUNE_THROW(Dune::NotImplemented,
178 "Optimized LocalOperator with matrix coefficients not implemented.");
185 static constexpr int degree = 0;
zero-order operator , or
Definition: ZeroOrderTestvecTrialvec.hpp:25
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: ZeroOrderTestvecTrialvec.hpp:19