5 #include <amdis/GridFunctionOperator.hpp> 6 #include <amdis/Output.hpp> 7 #include <amdis/common/StaticSize.hpp> 8 #include <amdis/common/ValueCategory.hpp> 24 template <
class LC,
class Gr
idFct>
32 using expr_value_type =
typename GridFct::Range;
33 static_assert( static_size_v<expr_value_type> == 1 || (static_num_rows_v<expr_value_type> == dow && static_num_cols_v<expr_value_type> == dow),
34 "Expression must be of scalar or matrix type." );
41 template <
class CG,
class RN,
class CN,
class Mat>
42 void getElementMatrix(CG
const& contextGeo, RN
const& rowNode, CN
const& colNode, Mat& elementMatrix)
44 static_assert(RN::isLeaf && CN::isLeaf,
45 "Operator can be applied to Leaf-Nodes only.");
47 const bool sameFE = std::is_same_v<FiniteElementType_t<RN>, FiniteElementType_t<CN>>;
48 const bool sameNode = rowNode.treeIndex() == colNode.treeIndex();
49 using Category = ValueCategory_t<typename GridFct::Range>;
51 auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
52 if (sameFE && sameNode)
53 getElementMatrixOptimized(contextGeo, quad, rowNode, colNode, elementMatrix, Category{});
55 getElementMatrixStandard(contextGeo, quad, rowNode, colNode, elementMatrix);
57 error_exit(
"Not implemented: currently only the implementation for equal fespaces available");
62 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
63 void getElementMatrixStandard(CG
const& contextGeo, QR
const& quad,
64 RN
const& rowNode, CN
const& colNode,
67 std::size_t size = rowNode.size();
69 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
70 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
71 std::vector<WorldVector> gradients;
73 for (
auto const& qp : quad) {
75 auto&& local = contextGeo.local(qp.position());
78 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
81 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
82 const auto exprValue = Super::coefficient(local);
85 auto const& shapeGradients = rowNode.localBasisJacobiansAt(local);
88 gradients.resize(shapeGradients.size());
90 for (std::size_t i = 0; i < gradients.size(); ++i)
91 jacobian.mv(shapeGradients[i][0], gradients[i]);
93 for (std::size_t i = 0; i < size; ++i) {
94 const auto local_i = rowNode.localIndex(i);
95 for (std::size_t j = 0; j < size; ++j) {
96 const auto local_j = colNode.localIndex(j);
97 elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
103 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
104 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
105 RN
const& node, CN
const& ,
108 std::size_t size = node.size();
110 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
111 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
112 std::vector<WorldVector> gradients;
114 for (
auto const& qp : quad) {
116 auto&& local = contextGeo.local(qp.position());
119 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
122 const auto factor = Super::coefficient(local) * contextGeo.integrationElement(qp.position())
126 auto const& shapeGradients = node.localBasisJacobiansAt(local);
129 gradients.resize(shapeGradients.size());
130 for (std::size_t i = 0; i < gradients.size(); ++i)
131 jacobian.mv(shapeGradients[i][0], gradients[i]);
133 for (std::size_t i = 0; i < size; ++i) {
134 const auto local_i = node.localIndex(i);
136 elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
138 for (std::size_t j = i+1; j < size; ++j) {
139 const auto local_j = node.localIndex(j);
140 const auto value = factor * (gradients[i] * gradients[j]);
142 elementMatrix[local_i][local_j] += value;
143 elementMatrix[local_j][local_i] += value;
149 template <
class CG,
class QR,
class RN,
class CN,
class Mat>
150 void getElementMatrixOptimized(CG
const& contextGeo, QR
const& quad,
151 RN
const& node, CN
const& ,
154 std::size_t size = node.size();
156 using RangeFieldType =
typename RN::LocalBasis::Traits::RangeFieldType;
157 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
158 std::vector<WorldVector> gradients;
160 for (
auto const& qp : quad) {
162 auto&& local = contextGeo.local(qp.position());
165 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
168 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
169 const auto exprValue = Super::coefficient(local);
172 auto const& shapeGradients = node.localBasisJacobiansAt(local);
175 gradients.resize(shapeGradients.size());
176 for (std::size_t i = 0; i < gradients.size(); ++i)
177 jacobian.mv(shapeGradients[i][0], gradients[i]);
179 for (std::size_t i = 0; i < size; ++i) {
180 const auto local_i = node.localIndex(i);
181 for (std::size_t j = 0; j < size; ++j) {
182 const auto local_j = node.localIndex(j);
183 elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
191 template <
class S,
class F,
class T,
int dow,
192 std::enable_if_t<Category::Scalar<S>,
int> = 0>
193 T eval(S
const& scalar, F factor,
194 Dune::FieldVector<T,dow>
const& grad_test,
195 Dune::FieldVector<T,dow>
const& grad_trial)
const 197 return (scalar * factor) * (grad_test * grad_trial);
200 template <
class M,
class F,
class T,
int dow,
201 std::enable_if_t<Category::Matrix<M>,
int> = 0>
202 T eval(M
const& mat, F factor,
203 Dune::FieldVector<T,dow>
const& grad_test,
204 Dune::FieldVector<T,dow>
const& grad_trial)
const 206 return factor * (grad_test * (mat * grad_trial));
211 template <
class Expr,
class... QuadratureArgs>
212 auto sot(Expr&& expr, QuadratureArgs&&... args)
Definition: SecondOrderGradTestGradTrial.hpp:19
auto makeOperator(Tag tag, Expr &&expr, QuadratureArgs &&... args)
Store tag and expression into a PreGridFunctionOperator to create a GridFunctionOperator.
Definition: GridFunctionOperator.hpp:220
The base-template for GridFunctionOperators.
Definition: GridFunctionOperator.hpp:242
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
auto sot(Expr &&expr, QuadratureArgs &&... args)
Create a second-order term.
Definition: SecondOrderGradTestGradTrial.hpp:212
Wrapper class for element and geometry.
Definition: ContextGeometry.hpp:43
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39