Fix an issue in SecondOrderDivTestvecDivTrialvec.hpp that made the...

Fix an issue in SecondOrderDivTestvecDivTrialvec.hpp that made the computations of the operator wrong.

amdis/localoperators/SecondOrderDivTestvecDivTrialvec.hpp

@@ -89,10 +89,12 @@ namespace AMDiS
 
         for (std::size_t i = 0; i < rowSize; ++i) {
           for (std::size_t j = 0; j < colSize; ++j) {
-            for (std::size_t k = 0; k < rowNode.degree(); ++k) {
+              for (std::size_t k = 0; k < CG::dow; ++k) {
+                  for (std::size_t l = 0; l < CG::dow; ++l) {
                       const auto local_ki = rowNode.child(k).localIndex(i);
-              const auto local_kj = colNode.child(k).localIndex(j);
-              elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][k];
+                      const auto local_kj = colNode.child(l).localIndex(j);
+                      elementMatrix[local_ki][local_kj] += factor * rowGradients[i][k] * colGradients[j][l];
+                  }
               }
           }
         }
@@ -130,16 +132,18 @@ namespace AMDiS
         for (std::size_t i = 0; i < gradients.size(); ++i)
           jacobian.mv(shapeGradients[i][0], gradients[i]);
 
-        for (std::size_t k = 0; k < node.degree(); ++k) {
+        for (std::size_t k = 0; k < CG::dow; ++k) {
+          for (std::size_t k = 0; k < CG::dow; ++k) {
             for (std::size_t i = 0; i < size; ++i) {
               const auto local_ki = node.child(k).localIndex(i);
               const auto value = factor * gradients[i][k];
               elementMatrix[local_ki][local_ki] += value * gradients[i][k];
 
-            for (std::size_t j = i+1; j < size; ++j) {
-              const auto local_kj = node.child(k).localIndex(j);
-              elementMatrix[local_ki][local_kj] += value * gradients[j][k];
-              elementMatrix[local_kj][local_ki] += value * gradients[j][k];
+              for (std::size_t j = i + 1; j < size; ++j) {
+                const auto local_kj = node.child(l).localIndex(j);
+                elementMatrix[local_ki][local_kj] += value * gradients[j][l];
+                elementMatrix[local_kj][local_ki] += value * gradients[j][l];
+              }
             }
           }
         }