5 #include <amdis/GridFunctionOperator.hpp> 6 #include <amdis/Output.hpp> 7 #include <amdis/common/StaticSize.hpp> 23 template <
class LC,
class Gr
idFct>
30 static_assert( static_size_v<typename GridFct::Range> == 1,
"Expression must be of scalar type." );
37 template <
class CG,
class RN,
class CN,
class Mat>
38 void getElementMatrix(CG
const& contextGeo, RN
const& rowNode, CN
const& colNode, Mat& elementMatrix)
40 static_assert(RN::isPower && CN::isPower,
41 "Operator can be applied to Power-Nodes only.");
43 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
44 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
46 auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
47 if (sameFE && sameNode)
48 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix);
50 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
54 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
55 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
56 RN
const& rowNode, CN
const& colNode,
59 static const std::size_t CHILDREN = RN::CHILDREN;
60 static_assert( RN::CHILDREN == CN::CHILDREN,
"" );
62 std::size_t rowSize = rowNode.child(0).size();
63 std::size_t colSize = colNode.child(0).size();
65 using RowFieldType =
typename RN::ChildType::LocalBasis::Traits::RangeFieldType;
66 using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
67 std::vector<RowWorldVector> rowGradients;
69 using ColFieldType =
typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
70 using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
71 std::vector<ColWorldVector> colGradients;
73 for (
auto const& qp : quad) {
75 decltype(
auto) local = contextGeo.local(qp.position());
78 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
81 const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
84 auto const& rowShapeGradients = rowNode.child(0).localBasisJacobiansAt(local);
85 auto const& colShapeGradients = colNode.child(0).localBasisJacobiansAt(local);
88 rowGradients.resize(rowShapeGradients.size());
90 for (std::size_t i = 0; i < rowGradients.size(); ++i)
91 jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
93 colGradients.resize(colShapeGradients.size());
95 for (std::size_t i = 0; i < colGradients.size(); ++i)
96 jacobian.mv(colShapeGradients[i][0], colGradients[i]);
98 for (std::size_t i = 0; i < rowSize; ++i) {
99 for (std::size_t j = 0; j < colSize; ++j) {
100 for (std::size_t k = 0; k < CHILDREN; ++k) {
101 const auto local_ki = rowNode.child(k).localIndex(i);
102 const auto local_kj = colNode.child(k).localIndex(j);
103 elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][k];
111 template <
class CG,
class QR,
class Node,
class Mat>
112 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
113 Node
const& node, Node
const& ,
116 static const std::size_t CHILDREN = Node::CHILDREN;
118 std::size_t size = node.child(0).size();
120 using RangeFieldType =
typename Node::ChildType::LocalBasis::Traits::RangeFieldType;
121 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
122 std::vector<WorldVector> gradients;
124 for (
auto const& qp : quad) {
126 auto&& local = contextGeo.local(qp.position());
129 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
132 const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
135 auto const& shapeGradients = node.child(0).localBasisJacobiansAt(local);
138 gradients.resize(shapeGradients.size());
140 for (std::size_t i = 0; i < gradients.size(); ++i)
141 jacobian.mv(shapeGradients[i][0], gradients[i]);
143 for (std::size_t k = 0; k < CHILDREN; ++k) {
144 for (std::size_t i = 0; i < size; ++i) {
145 const auto local_ki = node.child(k).localIndex(i);
146 const auto value = factor * gradients[i][k];
147 elementMatrix[local_ki][local_ki] += value * gradients[i][k];
149 for (std::size_t j = i+1; j < size; ++j) {
150 const auto local_kj = node.child(k).localIndex(j);
151 elementMatrix[local_ki][local_kj] += value * gradients[j][k];
152 elementMatrix[local_kj][local_ki] += value * gradients[j][k];
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
Definition: SecondOrderDivTestvecDivTrialvec.hpp:18
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39