AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
SecondOrderDivTestvecDivTrialvec.hpp
1#pragma once
2
3#include <type_traits>
4
5#include <amdis/GridFunctionOperator.hpp>
6#include <amdis/common/StaticSize.hpp>
7#include <amdis/typetree/FiniteElementType.hpp>
8
9namespace AMDiS
10{
16 namespace tag
17 {
19 }
20
21
24 {
25 public:
27
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
31 {
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.");
36
37 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
38 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
39
40 if (sameFE && sameNode)
41 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
42 else
43 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, localFct, elementMatrix);
44 }
45
46 protected:
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
51 {
52 assert( rowNode.degree() == colNode.degree() );
53
54 std::size_t rowSize = rowNode.child(0).size();
55 std::size_t colSize = colNode.child(0).size();
56
57 using RowFieldType = typename RN::ChildType::LocalBasis::Traits::RangeFieldType;
58 using RowWorldVector = FieldVector<RowFieldType,CG::dow>;
59 std::vector<RowWorldVector> rowGradients;
60
61 using ColFieldType = typename CN::ChildType::LocalBasis::Traits::RangeFieldType;
62 using ColWorldVector = FieldVector<ColFieldType,CG::dow>;
63 std::vector<ColWorldVector> colGradients;
64
65 for (auto const& qp : quad) {
66 // Position of the current quadrature point in the reference element
67 decltype(auto) local = contextGeo.coordinateInElement(qp.position());
68
69 // The transposed inverse Jacobian of the map from the reference element to the element
70 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
71
72 // The multiplicative factor in the integral transformation formula
73 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
74
75 // The gradients of the shape functions on the reference element
76 auto const& rowShapeGradients = rowNode.child(0).localBasisJacobiansAt(local);
77 auto const& colShapeGradients = colNode.child(0).localBasisJacobiansAt(local);
78
79 // Compute the shape function gradients on the real element
80 rowGradients.resize(rowShapeGradients.size());
81
82 for (std::size_t i = 0; i < rowGradients.size(); ++i)
83 jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
84
85 colGradients.resize(colShapeGradients.size());
86
87 for (std::size_t i = 0; i < colGradients.size(); ++i)
88 jacobian.mv(colShapeGradients[i][0], colGradients[i]);
89
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];
97 }
98 }
99 }
100 }
101 }
102 }
103
104
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& /*colNode*/,
108 LocalFct const& localFct, Mat& elementMatrix) const
109 {
110 std::size_t size = node.child(0).size();
111
112 using RangeFieldType = typename Node::ChildType::LocalBasis::Traits::RangeFieldType;
113 using WorldVector = FieldVector<RangeFieldType,CG::dow>;
114 std::vector<WorldVector> gradients;
115
116 for (auto const& qp : quad) {
117 // Position of the current quadrature point in the reference element
118 auto&& local = contextGeo.coordinateInElement(qp.position());
119
120 // The transposed inverse Jacobian of the map from the reference element to the element
121 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
122
123 // The multiplicative factor in the integral transformation formula
124 const auto factor = localFct(local) * contextGeo.integrationElement(qp.position()) * qp.weight();
125
126 // The gradients of the shape functions on the reference element
127 auto const& shapeGradients = node.child(0).localBasisJacobiansAt(local);
128
129 // Compute the shape function gradients on the real element
130 gradients.resize(shapeGradients.size());
131
132 for (std::size_t i = 0; i < gradients.size(); ++i)
133 jacobian.mv(shapeGradients[i][0], gradients[i]);
134
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];
141
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];
146 }
147 }
148 }
149 }
150 }
151 }
152 };
153
154 template <class LC>
155 struct GridFunctionOperatorRegistry<tag::divtestvec_divtrialvec, LC>
156 {
157 static constexpr int degree = 2;
159 };
160
163} // end namespace AMDiS
second-order operator
Definition: SecondOrderDivTestvecDivTrialvec.hpp:24
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: SecondOrderDivTestvecDivTrialvec.hpp:18