discretization: fix evaluation of stencil grantient weights

src/containers/piercedKernelIterator.hpp

@@ -130,17 +130,15 @@ friend class PiercedStorageIterator;
     /**
     * Two-way comparison.
     */
-    template<typename other_id_t>
-    bool operator==(const PiercedKernelIterator<other_id_t>& rhs) const
+    bool operator==(const PiercedKernelIterator &rhs) const
     {
         return (m_kernel == rhs.m_kernel) && (m_pos == rhs.m_pos);
     }
 
     /**
     * Two-way comparison.
     */
-    template<typename other_id_t>
-    bool operator!=(const PiercedKernelIterator<other_id_t>& rhs) const
+    bool operator!=(const PiercedKernelIterator &rhs) const
     {
         return (m_kernel != rhs.m_kernel) || (m_pos != rhs.m_pos);
     }

src/containers/piercedStorageIterator.hpp

@@ -160,19 +160,15 @@ friend class PiercedStorage;
     /**
     * Two-way comparison.
     */
-    template<typename other_value_t, typename other_id_t = long,
-         typename other_value_no_cv_t = typename std::remove_cv<other_value_t>::type>
-    bool operator==(const PiercedStorageIterator<other_value_t, other_id_t, other_value_no_cv_t>& rhs) const
+    bool operator==(const PiercedStorageIterator &rhs) const
     {
         return (PiercedKernelIterator<id_t>::operator==(rhs) && m_storage == rhs.m_storage);
     }
 
     /**
     * Two-way comparison.
     */
-    template<typename other_value_t, typename other_id_t = long,
-         typename other_value_no_cv_t = typename std::remove_cv<other_value_t>::type>
-    bool operator!=(const PiercedStorageIterator<other_value_t, other_id_t, other_value_no_cv_t>& rhs) const
+    bool operator!=(const PiercedStorageIterator &rhs) const
     {
         return (PiercedKernelIterator<id_t>::operator!=(rhs) || m_storage != rhs.m_storage);
     }

src/discretization/reconstruction.cpp

@@ -2673,51 +2673,34 @@ void ReconstructionKernel::computeGradientLimitedWeights(uint8_t degree, const s
     int nCoeffs = ReconstructionPolynomial::getCoefficientCount(degree, dimensions);
 
     BITPIT_CREATE_WORKSPACE(dcsi, double, nCoeffs, MAX_STACK_WORKSPACE_SIZE);
-    BITPIT_CREATE_WORKSPACE(directionalStencilWeights, double, nEquations, MAX_STACK_WORKSPACE_SIZE);
 
     const double *polynomialWeights = getPolynomialWeights();
 
-    // Evalaute weights in x-direction
-    ReconstructionPolynomial::evalPointBasisDerivatives(degree, dimensions, origin, point, {{1., 0., 0.}}, dcsi);
-    if (limiters) {
-        applyLimiter(degree, limiters, dcsi);
-    }
-
-    cblas_dgemv(CBLAS_ORDER::CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans,
-                nEquations, nCoeffs, 1., polynomialWeights, nEquations,
-                dcsi, 1, 0, directionalStencilWeights, 1);
-
-    for (int i = 0; i < nEquations; ++i) {
-        gradientWeights[i][0] = directionalStencilWeights[i];
-    }
-
-    // Evalaute weights in y-direction
-    ReconstructionPolynomial::evalPointBasisDerivatives(degree, getDimensions(), origin, point, {{0., 1., 0.}}, dcsi);
-    if (limiters) {
-        applyLimiter(degree, limiters, dcsi);
-    }
+    std::array<double, 3> direction = {{0., 0., 0.}};
+    for (int d = 0; d < dimensions; ++d) {
+        direction[d] = 1.;
 
-    cblas_dgemv(CBLAS_ORDER::CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans,
-                nEquations, nCoeffs, 1., polynomialWeights, nEquations,
-                dcsi, 1, 0, directionalStencilWeights, 1);
-
-    for (int i = 0; i < nEquations; ++i) {
-        gradientWeights[i][1] = directionalStencilWeights[i];
-    }
-
-    // Evalaute weights in z-direction
-    if (dimensions == 3) {
-        ReconstructionPolynomial::evalPointBasisDerivatives(degree, getDimensions(), origin, point, {{0., 0., 1.}}, dcsi);
+        ReconstructionPolynomial::evalPointBasisDerivatives(degree, dimensions, origin, point, direction, dcsi);
         if (limiters) {
             applyLimiter(degree, limiters, dcsi);
         }
 
+        // Weights are stored in contiguous three-dimensional arrays, this
+        // means we can access the weights of the current dimensions using
+        // the dimension as offset and a stride of three elements
         cblas_dgemv(CBLAS_ORDER::CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans,
                     nEquations, nCoeffs, 1., polynomialWeights, nEquations,
-                    dcsi, 1, 0, directionalStencilWeights, 1);
+                    dcsi, 1, 0, gradientWeights->data() + d, 3);
+
+        direction[d] = 0.;
+    }
 
+    // Explicitly zero unused components
+    if (dimensions != ReconstructionPolynomial::MAX_DIMENSIONS) {
         for (int i = 0; i < nEquations; ++i) {
-            gradientWeights[i][2] = directionalStencilWeights[i];
+            for (int d = dimensions; d < ReconstructionPolynomial::MAX_DIMENSIONS; ++d) {
+                gradientWeights[i][d] = 0.;
+            }
         }
     }
 }