5 #include <amdis/GridFunctionOperator.hpp> 6 #include <amdis/common/StaticSize.hpp> 7 #include <amdis/typetree/FiniteElementType.hpp> 28 template <
class CG,
class RN,
class CN,
class Quad,
class LocalFct,
class Mat>
29 void assemble(CG
const& contextGeo, RN
const& rowNode, CN
const& colNode,
30 Quad
const& quad, LocalFct
const& localFct, Mat& elementMatrix)
const 32 static_assert(static_size_v<typename LocalFct::Range> == 1,
33 "Expression must be of scalar type.");
34 static_assert(RN::isPower && CN::isPower,
35 "Operator can be applied to Power-Nodes only.");
37 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
38 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
40 if (sameFE && sameNode)
41 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
43 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
47 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
48 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
49 RN
const& rowNode, CN
const& colNode,
50 LocalFct
const& localFct, Mat& elementMatrix)
const 52 static const std::size_t CHILDREN = RN::CHILDREN;
53 static_assert( RN::CHILDREN == CN::CHILDREN,
"" );
55 std::size_t rowSize = rowNode.child(0).size();
56 std::size_t colSize = colNode.child(0).size();
58 using RowFieldType =
typename RN::ChildType::LocalBasis::Traits::RangeFieldType;
59 using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
60 std::vector<RowWorldVector> rowGradients;
62 using ColFieldType =
typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
63 using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
64 std::vector<ColWorldVector> colGradients;
66 for (
auto const& qp : quad) {
68 decltype(
auto) local = contextGeo.local(qp.position());
71 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
74 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
77 auto const& rowShapeGradients = rowNode.child(0).localBasisJacobiansAt(local);
78 auto const& colShapeGradients = colNode.child(0).localBasisJacobiansAt(local);
81 rowGradients.resize(rowShapeGradients.size());
83 for (std::size_t i = 0; i < rowGradients.size(); ++i)
84 jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
86 colGradients.resize(colShapeGradients.size());
88 for (std::size_t i = 0; i < colGradients.size(); ++i)
89 jacobian.mv(colShapeGradients[i][0], colGradients[i]);
91 for (std::size_t i = 0; i < rowSize; ++i) {
92 for (std::size_t j = 0; j < colSize; ++j) {
93 for (std::size_t k = 0; k < CHILDREN; ++k) {
94 const auto local_ki = rowNode.child(k).localIndex(i);
95 const auto local_kj = colNode.child(k).localIndex(j);
96 elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][k];
104 template <
class CG,
class QR,
class Node,
class LocalFct,
class Mat>
105 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
106 Node
const& node, Node
const& ,
107 LocalFct
const& localFct, Mat& elementMatrix)
const 109 static const std::size_t CHILDREN = Node::CHILDREN;
111 std::size_t size = node.child(0).size();
113 using RangeFieldType =
typename Node::ChildType::LocalBasis::Traits::RangeFieldType;
114 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
115 std::vector<WorldVector> gradients;
117 for (
auto const& qp : quad) {
119 auto&& local = contextGeo.local(qp.position());
122 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
125 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
128 auto const& shapeGradients = node.child(0).localBasisJacobiansAt(local);
131 gradients.resize(shapeGradients.size());
133 for (std::size_t i = 0; i < gradients.size(); ++i)
134 jacobian.mv(shapeGradients[i][0], gradients[i]);
136 for (std::size_t k = 0; k < CHILDREN; ++k) {
137 for (std::size_t i = 0; i < size; ++i) {
138 const auto local_ki = node.child(k).localIndex(i);
139 const auto value = factor * gradients[i][k];
140 elementMatrix[local_ki][local_ki] += value * gradients[i][k];
142 for (std::size_t j = i+1; j < size; ++j) {
143 const auto local_kj = node.child(k).localIndex(j);
144 elementMatrix[local_ki][local_kj] += value * gradients[j][k];
145 elementMatrix[local_kj][local_ki] += value * gradients[j][k];
156 static constexpr
int degree = 2;
Definition: AdaptBase.hpp:6
Definition: SecondOrderDivTestvecDivTrialvec.hpp:18
second-order operator
Definition: SecondOrderDivTestvecDivTrialvec.hpp:23
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:204