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 static_assert(RN::CHILDREN == CN::CHILDREN);
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 static const std::size_t CHILDREN = RN::CHILDREN;
62 std::size_t rowSize = rowNode.child(0).size();
63 std::size_t colSize = colNode.child(0).size();
65 for (
auto const& qp : quad) {
67 auto&& local = contextGeo.local(qp.position());
70 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
72 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
73 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
75 for (std::size_t i = 0; i < rowSize; ++i) {
76 const auto value = factor * rowShapeValues[i];
78 for (std::size_t j = 0; j < colSize; ++j) {
79 const auto value0 = value * colShapeValues[j];
81 for (std::size_t k = 0; k < CHILDREN; ++k) {
82 const auto local_ki = rowNode.child(k).localIndex(i);
83 const auto local_kj = colNode.child(k).localIndex(j);
85 elementMatrix[local_ki][local_kj] += value0;
93 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
94 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
95 RN
const& node, CN
const& ,
96 LocalFct
const& localFct, Mat& elementMatrix,
tag::scalar)
const 98 static const std::size_t CHILDREN = RN::CHILDREN;
99 std::size_t size = node.child(0).size();
101 for (
auto const& qp : quad) {
103 auto&& local = contextGeo.local(qp.position());
106 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
108 auto const& shapeValues = node.child(0).localBasisValuesAt(local);
110 for (std::size_t i = 0; i < size; ++i) {
111 const auto value = factor * shapeValues[i];
113 for (std::size_t k = 0; k < CHILDREN; ++k) {
114 const auto local_ki = node.child(k).localIndex(i);
115 elementMatrix[local_ki][local_ki] += value * shapeValues[i];
118 for (std::size_t j = i+1; j < size; ++j) {
119 const auto value0 = value * shapeValues[j];
121 for (std::size_t k = 0; k < CHILDREN; ++k) {
122 const auto local_ki = node.child(k).localIndex(i);
123 const auto local_kj = node.child(k).localIndex(j);
125 elementMatrix[local_ki][local_kj] += value0;
126 elementMatrix[local_kj][local_ki] += value0;
133 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
134 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
135 RN
const& rowNode, CN
const& colNode,
136 LocalFct
const& localFct, Mat& elementMatrix,
tag::matrix)
const 138 static const std::size_t CHILDREN = RN::CHILDREN;
139 using expr_value_type =
typename LocalFct::Range;
140 static_assert(static_num_rows_v<expr_value_type> == CHILDREN &&
141 static_num_cols_v<expr_value_type> == CHILDREN);
143 std::size_t rowSize = rowNode.child(0).size();
144 std::size_t colSize = colNode.child(0).size();
146 for (
auto const& qp : quad) {
148 auto&& local = contextGeo.local(qp.position());
151 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
152 const auto exprValue = localFct(local);
154 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
155 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
157 for (std::size_t i = 0; i < rowSize; ++i) {
158 const auto value0 = factor * rowShapeValues[i];
160 for (std::size_t j = 0; j < colSize; ++j) {
161 const auto value = value0 * colShapeValues[j];
162 const auto mat = exprValue * value;
164 for (std::size_t k0 = 0; k0 < CHILDREN; ++k0) {
165 const auto local_ki = rowNode.child(k0).localIndex(i);
166 for (std::size_t k1 = 0; k1 < CHILDREN; ++k1) {
167 const auto local_kj = colNode.child(k1).localIndex(j);
169 elementMatrix[local_ki][local_kj] += mat[k0][k1];
177 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
178 void getElementMatrixOptimized(CG
const&, QR
const&, RN
const&, CN
const&,
181 DUNE_THROW(Dune::NotImplemented,
182 "Optimized LocalOperator with matrix coefficients not implemented.");
189 static constexpr
int degree = 0;
Definition: AdaptBase.hpp:6
zero-order operator , or
Definition: ZeroOrderTestvecTrialvec.hpp:24
Definition: ZeroOrderTestvecTrialvec.hpp:19
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:204