6#include <amdis/GridFunctionOperator.hpp>
7#include <amdis/common/StaticSize.hpp>
20 bool constViscosity =
false;
32 : constViscosity_(t.constViscosity)
35 template <
class CG,
class Node,
class Quad,
class LocalFct,
class Mat>
36 void assemble(CG
const& contextGeo, Node
const& tree, Node
const& ,
37 Quad
const& quad, LocalFct
const& localFct, Mat& elementMatrix)
const
39 static_assert(static_size_v<typename LocalFct::Range> == 1,
40 "Viscosity must be of scalar type.");
42 using VelocityNode = Dune::TypeTree::Child<Node,0>;
43 using VelocityCompNode = Dune::TypeTree::Child<Node,0,0>;
44 using PressureNode = Dune::TypeTree::Child<Node,1>;
46 static_assert(VelocityNode::isPower && PressureNode::isLeaf,
"");
48 using namespace Dune::Indices;
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.coordinateInElement(qp.position());
62 const auto jacobian = contextGeo.elementGeometry().jacobianInverseTransposed(local);
65 const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
66 const auto vel_factor = localFct(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;
133 static constexpr int degree = 2;
Stokes operator .
Definition: StokesOperator.hpp:27
Registry to specify a tag for each implementation type.
Definition: GridFunctionOperator.hpp:216
Definition: StokesOperator.hpp:19