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." );
35 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
36 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
38 if (sameFE && sameNode)
39 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
41 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
44 template <
class CG,
class Node,
class Quad,
class LocalFct,
class Vec>
45 void assemble(CG
const& contextGeo, Node
const& node,
46 Quad
const& quad, LocalFct
const& localFct, Vec& elementVector)
const
48 static_assert(static_size_v<typename LocalFct::Range> == 1,
49 "Expression must be of scalar type." );
50 static_assert(Node::isLeaf,
51 "Operator can be applied to Leaf-Nodes only");
53 std::size_t size = node.size();
55 for (
auto const& qp : quad) {
57 auto&& local = contextGeo.coordinateInElement(qp.position());
60 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
62 auto const& shapeValues = node.localBasisValuesAt(local);
63 for (std::size_t i = 0; i < size; ++i) {
64 const auto local_i = node.localIndex(i);
65 elementVector[local_i] += factor * shapeValues[i];
72 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
73 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
74 RN
const& rowNode, CN
const& colNode,
75 LocalFct
const& localFct, Mat& elementMatrix)
const
77 static_assert(RN::isLeaf && CN::isLeaf,
78 "Operator can be applied to Leaf-Nodes only.");
80 std::size_t rowSize = rowNode.size();
81 std::size_t colSize = colNode.size();
83 for (
auto const& qp : quad) {
85 auto&& local = contextGeo.coordinateInElement(qp.position());
88 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
90 auto const& rowShapeValues = rowNode.localBasisValuesAt(local);
91 auto const& colShapeValues = colNode.localBasisValuesAt(local);
93 for (std::size_t i = 0; i < rowSize; ++i) {
94 const auto local_i = rowNode.localIndex(i);
95 const auto value = factor * rowShapeValues[i];
97 for (std::size_t j = 0; j < colSize; ++j) {
98 const auto local_j = colNode.localIndex(j);
99 elementMatrix[local_i][local_j] += value * colShapeValues[j];
106 template <
class CG,
class QR,
class RN,
class CN,
class LocalFct,
class Mat>
107 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
108 RN
const& node, CN
const& ,
109 LocalFct
const& localFct, Mat& elementMatrix)
const
111 static_assert(RN::isLeaf && CN::isLeaf,
112 "Operator can be applied to Leaf-Nodes only.");
114 std::size_t size = node.size();
116 for (
auto const& qp : quad) {
118 auto&& local = contextGeo.coordinateInElement(qp.position());
121 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
123 auto const& shapeValues = node.localBasisValuesAt(local);
125 for (std::size_t i = 0; i < size; ++i) {
126 const auto local_i = node.localIndex(i);
128 const auto value = factor * shapeValues[i];
129 elementMatrix[local_i][local_i] += value * shapeValues[i];
131 for (std::size_t j = i+1; j < size; ++j) {
132 const auto local_j = node.localIndex(j);
134 elementMatrix[local_i][local_j] += value * shapeValues[j];
135 elementMatrix[local_j][local_i] += value * shapeValues[j];
145 static constexpr int degree = 0;
151 template <
class Expr>
152 auto zot(Expr&& expr,
int quadOrder = -1)
zero-order operator
Definition: ZeroOrderTestTrial.hpp:24
auto makeOperator(Tag const &tag, Expr &&expr, int gridFctDeg=-1)
Definition: GridFunctionOperator.hpp:235
auto zot(Expr &&expr, int quadOrder=-1)
Create a zero-order term.
Definition: ZeroOrderTestTrial.hpp:152
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: ZeroOrderTestTrial.hpp:18