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 assert( rowNode.degree() == colNode.degree() );
54 std::size_t rowSize = rowNode.child(0).size();
55 std::size_t colSize = colNode.child(0).size();
57 using RowFieldType =
typename RN::ChildType::LocalBasis::Traits::RangeFieldType;
58 using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
59 std::vector<RowWorldVector> rowGradients;
61 using ColFieldType =
typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
62 using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
63 std::vector<ColWorldVector> colGradients;
65 for (
auto const& qp : quad) {
67 decltype(
auto) local = contextGeo.coordinateInElement(qp.position());
70 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
73 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
76 auto const& rowShapeGradients = rowNode.child(0).localBasisJacobiansAt(local);
77 auto const& colShapeGradients = colNode.child(0).localBasisJacobiansAt(local);
80 rowGradients.resize(rowShapeGradients.size());
82 for (std::size_t i = 0; i < rowGradients.size(); ++i)
83 jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
85 colGradients.resize(colShapeGradients.size());
87 for (std::size_t i = 0; i < colGradients.size(); ++i)
88 jacobian.mv(colShapeGradients[i][0], colGradients[i]);
90 for (std::size_t i = 0; i < rowSize; ++i) {
91 for (std::size_t j = 0; j < colSize; ++j) {
92 for (std::size_t k = 0; k < CG::dow; ++k) {
93 for (std::size_t l = 0; l < CG::dow; ++l) {
94 const auto local_ki = rowNode.child(k).localIndex(i);
95 const auto local_kj = colNode.child(l).localIndex(j);
96 elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][l];
105 template <
class CG,
class QR,
class Node,
class LocalFct,
class Mat>
106 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
107 Node
const& node, Node
const& ,
108 LocalFct
const& localFct, Mat& elementMatrix)
const
110 std::size_t size = node.child(0).size();
112 using RangeFieldType =
typename Node::ChildType::LocalBasis::Traits::RangeFieldType;
113 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
114 std::vector<WorldVector> gradients;
116 for (
auto const& qp : quad) {
118 auto&& local = contextGeo.coordinateInElement(qp.position());
121 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
124 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
127 auto const& shapeGradients = node.child(0).localBasisJacobiansAt(local);
130 gradients.resize(shapeGradients.size());
132 for (std::size_t i = 0; i < gradients.size(); ++i)
133 jacobian.mv(shapeGradients[i][0], gradients[i]);
135 for (std::size_t k = 0; k < CG::dow; ++k) {
136 for (std::size_t l = 0; l < CG::dow; ++l) {
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(l).localIndex(j);
144 elementMatrix[local_ki][local_kj] += value * gradients[j][l];
145 elementMatrix[local_kj][local_ki] += value * gradients[j][l];
157 static constexpr int degree = 2;
second-order operator
Definition: SecondOrderDivTestvecDivTrialvec.hpp:24
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: SecondOrderDivTestvecDivTrialvec.hpp:18