5 #include <amdis/GridFunctionOperator.hpp> 6 #include <amdis/common/StaticSize.hpp> 7 #include <amdis/common/ValueCategory.hpp> 23 template <
class LC,
class Gr
idFct>
31 using expr_value_type =
typename GridFct::Range;
32 static_assert( static_size_v<expr_value_type> == 1 || (static_num_rows_v<expr_value_type> == dow && static_num_cols_v<expr_value_type> == dow),
33 "Expression must be of scalar or matrix type." );
40 template <
class CG,
class RN,
class CN,
class Mat>
41 void getElementMatrix(CG
const& contextGeo, RN
const& rowNode, CN
const& colNode, Mat& elementMatrix)
43 static_assert(RN::isPower && CN::isPower,
44 "Operator can be applied to Power-Nodes only.");
46 static_assert( (RN::CHILDREN == CN::CHILDREN),
"" );
49 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
50 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
52 using Category = ValueCategory_t<expr_value_type>;
54 auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
55 if (sameFE && sameNode && std::is_same_v<Category,tag::scalar>)
56 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
58 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
63 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
64 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
65 RN
const& rowNode, CN
const& colNode,
68 static const std::size_t CHILDREN = RN::CHILDREN;
70 std::size_t rowSize = rowNode.child(0).size();
71 std::size_t colSize = colNode.child(0).size();
73 for (
auto const& qp : quad) {
75 auto&& local = contextGeo.local(qp.position());
78 const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
80 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
81 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
83 for (std::size_t i = 0; i < rowSize; ++i) {
84 const auto value = factor * rowShapeValues[i];
86 for (std::size_t j = 0; j < colSize; ++j) {
87 const auto value0 = value * colShapeValues[j];
89 for (std::size_t k = 0; k < CHILDREN; ++k) {
90 const auto local_ki = rowNode.child(k).localIndex(i);
91 const auto local_kj = colNode.child(k).localIndex(j);
93 elementMatrix[local_ki][local_kj] += value0;
101 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
102 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
103 RN
const& node, CN
const& ,
106 static const std::size_t CHILDREN = RN::CHILDREN;
107 std::size_t size = node.child(0).size();
109 for (
auto const& qp : quad) {
111 auto&& local = contextGeo.local(qp.position());
114 const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
116 auto const& shapeValues = node.child(0).localBasisValuesAt(local);
118 for (std::size_t i = 0; i < size; ++i) {
119 const auto value = factor * shapeValues[i];
121 for (std::size_t k = 0; k < CHILDREN; ++k) {
122 const auto local_ki = node.child(k).localIndex(i);
123 elementMatrix[local_ki][local_ki] += value * shapeValues[i];
126 for (std::size_t j = i+1; j < size; ++j) {
127 const auto value0 = value * shapeValues[j];
129 for (std::size_t k = 0; k < CHILDREN; ++k) {
130 const auto local_ki = node.child(k).localIndex(i);
131 const auto local_kj = node.child(k).localIndex(j);
133 elementMatrix[local_ki][local_kj] += value0;
134 elementMatrix[local_kj][local_ki] += value0;
141 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
142 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
143 RN
const& rowNode, CN
const& colNode,
146 static const std::size_t CHILDREN = RN::CHILDREN;
147 static_assert(static_num_rows_v<expr_value_type> == CHILDREN && static_num_cols_v<expr_value_type> == CHILDREN,
"" );
149 std::size_t rowSize = rowNode.child(0).size();
150 std::size_t colSize = colNode.child(0).size();
152 for (
auto const& qp : quad) {
154 auto&& local = contextGeo.local(qp.position());
157 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
158 const auto exprValue = Super::coefficient(local);
160 auto const& rowShapeValues = rowNode.child(0).localBasisValuesAt(local);
161 auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);
163 for (std::size_t i = 0; i < rowSize; ++i) {
164 const auto value0 = factor * rowShapeValues[i];
166 for (std::size_t j = 0; j < colSize; ++j) {
167 const auto value = value0 * colShapeValues[j];
168 const auto mat = exprValue * value;
170 for (std::size_t k0 = 0; k0 < CHILDREN; ++k0) {
171 const auto local_ki = rowNode.child(k0).localIndex(i);
172 for (std::size_t k1 = 0; k1 < CHILDREN; ++k1) {
173 const auto local_kj = colNode.child(k1).localIndex(j);
175 elementMatrix[local_ki][local_kj] += mat[k0][k1];
183 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
184 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
185 RN
const& node, CN
const& ,
188 assert(
false &&
"Not implemented.");
The base-template for GridFunctionOperators.
Definition: GridFunctionOperator.hpp:242
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Wrapper class for element and geometry.
Definition: ContextGeometry.hpp:43
Definition: ZeroOrderTestvecTrialvec.hpp:18
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39