6 #include <amdis/GridFunctionOperator.hpp> 7 #include <amdis/common/StaticSize.hpp> 18 struct stokes {
bool constViscosity =
false; };
23 template <
class LC,
class ViscosityExpr>
30 static_assert( static_size_v<typename ViscosityExpr::Range> == 1,
"Viscosity must be of scalar type." );
35 , constViscosity_(t.constViscosity)
38 template <
class CG,
class Node,
class Mat>
39 void getElementMatrix(CG
const& contextGeo, Node
const& tree, Node
const& , Mat& elementMatrix)
41 using VelocityNode =
typename Node::template Child<0>::type;
42 using VelocityCompNode =
typename VelocityNode::template Child<0>::type;
43 using PressureNode =
typename Node::template Child<1>::type;
44 static_assert(VelocityNode::isPower && PressureNode::isLeaf,
"");
46 using namespace Dune::Indices;
48 auto const& quad = this->getQuadratureRule(contextGeo.type(), tree, tree);
50 std::size_t numVelocityLocalFE = tree.child(_0,0).size();
51 std::size_t numPressureLocalFE = tree.child(_1).size();
53 using RangeFieldType =
typename VelocityNode::ChildType::LocalBasis::Traits::RangeFieldType;
54 using WorldVector = Dune::FieldVector<RangeFieldType,CG::dow>;
55 std::vector<WorldVector> gradients;
57 for (
auto const& qp : quad) {
59 auto&& local = contextGeo.local(qp.position());
62 const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
65 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
66 const auto vel_factor = Super::coefficient(local) * factor;
68 auto const& shapeGradients = tree.child(_0,0).localBasisJacobiansAt(local);
71 gradients.resize(shapeGradients.size());
72 for (std::size_t i = 0; i < gradients.size(); ++i)
73 jacobian.mv(shapeGradients[i][0], gradients[i]);
75 auto const& pressureValues = tree.child(_1).localBasisValuesAt(local);
78 for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
79 const auto value_ii = vel_factor * (gradients[i] * gradients[i]);
80 for (std::size_t k = 0; k < CG::dow; ++k) {
81 const auto local_ki = tree.child(_0,k).localIndex(i);
82 elementMatrix[local_ki][local_ki] += value_ii;
85 for (std::size_t j = i+1; j < numVelocityLocalFE; ++j) {
86 const auto value = vel_factor * (gradients[i] * gradients[j]);
87 for (std::size_t k = 0; k < CG::dow; ++k) {
88 const auto local_ki = tree.child(_0,k).localIndex(i);
89 const auto local_kj = tree.child(_0,k).localIndex(j);
90 elementMatrix[local_ki][local_kj] += value;
91 elementMatrix[local_kj][local_ki] += value;
96 if (!constViscosity_) {
98 for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
99 for (std::size_t kj = 0; kj < CG::dow; ++kj) {
100 const auto value_i_kj = vel_factor * gradients[i][kj];
101 for (std::size_t j = 0; j < numVelocityLocalFE; ++j) {
102 for (std::size_t ki = 0; ki < CG::dow; ++ki) {
103 const auto local_ki = tree.child(_0,ki).localIndex(i);
104 const auto local_kj = tree.child(_0,kj).localIndex(j);
105 elementMatrix[local_ki][local_kj] += value_i_kj * gradients[j][ki];
113 for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
114 for (std::size_t j = 0; j < numPressureLocalFE; ++j) {
115 const auto value = factor * pressureValues[j];
116 for (std::size_t k = 0; k < CG::dow; ++k) {
117 const auto local_v = tree.child(_0,k).localIndex(i);
118 const auto local_p = tree.child(_1).localIndex(j);
120 elementMatrix[local_v][local_p] += gradients[i][k] * value;
121 elementMatrix[local_p][local_v] += gradients[i][k] * value;
130 bool constViscosity_;
The base-template for GridFunctionOperators.
Definition: GridFunctionOperator.hpp:242
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Definition: StokesOperator.hpp:18
The main implementation of an CRTP-base class for operators using a grid-function coefficient to be u...
Definition: GridFunctionOperator.hpp:39