added stress tensor part for non-constant viscosity

src/amdis/localoperators/StokesOperator.hpp

@@ -31,8 +31,9 @@ namespace AMDiS
     static_assert( Size_v<typename ViscosityExpr::Range> == 1, "Viscosity must be of scalar type." );
 
   public:
-    GridFunctionOperator(tag::stokes, ViscosityExpr const& expr)
+    GridFunctionOperator(tag::stokes, ViscosityExpr const& expr, bool constViscosity = false)
       : Super(expr, 1)
+      , constViscosity_(constViscosity)
     {}
 
     template <class CG, class Node, class Mat>
@@ -100,6 +101,22 @@ namespace AMDiS
           }
         }
 
+        if (!constViscosity_) {
+          // <viscosity * d_i u_j, d_j v_i>
+          for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
+          for (std::size_t kj = 0; kj < CG::dow; ++kj) {
+            const auto value_i_kj = vel_factor * gradients[i][kj];
+            for (std::size_t j = 0; j < numVelocityLocalFE; ++j) {
+            for (std::size_t ki = 0; ki < CG::dow; ++ki) {
+              const auto local_ki = tree.child(_0,ki).localIndex(i);
+              const auto local_kj = tree.child(_0,kj).localIndex(j);
+              elementMatrix[local_ki][local_kj] += value_i_kj * gradients[j][ki];
+            }
+            }
+          }
+          }
+        }
+
         // <p, div(v)> + <div(u), q>
         for (std::size_t i = 0; i < numVelocityLocalFE; ++i) {
           for (std::size_t j = 0; j < numPressureLocalFE; ++j) {
@@ -116,6 +133,9 @@ namespace AMDiS
       }
 
     } // end calculateElementMatrix
+
+  private:
+    bool constViscosity_;
   };
 
   /** @} **/