ViennaTools · FilipovicLado · Jan 15, 2026 · Dec 6, 2025 · Dec 7, 2025 · Dec 9, 2025
diff --git a/.gitignore b/.gitignore
@@ -23,4 +23,5 @@ _generate/
 *env*
 .mypy_cache/
 .eggs/
-test.ipynb
+test.ipynb
+testall.sh
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -112,7 +112,7 @@ include(cmake/vtk.cmake)
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 1.9.0
+  VERSION 1.9.3
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   OPTIONS "VIENNACORE_FORMAT_EXCLUDE build/"
   EXCLUDE_FROM_ALL ${VIENNALS_BUILD_PYTHON})
@@ -125,7 +125,7 @@ CPMAddPackage(
 
 CPMFindPackage(
   NAME ViennaHRLE
-  VERSION 0.7.0
+  VERSION 0.8.0
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaHRLE"
   EXCLUDE_FROM_ALL ${VIENNALS_BUILD_PYTHON})
 

diff --git a/examples/SquareEtch/SquareEtch.cpp b/examples/SquareEtch/SquareEtch.cpp
@@ -200,8 +200,7 @@ int main() {
     // Analytical velocity fields and dissipation coefficients
     // can only be used with this spatial discretization scheme
     advectionKernel.setSpatialScheme(
-        // ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
-        ls::SpatialSchemeEnum::WENO_5TH_ORDER);
+        ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
   } else {
     // for numerical velocities, just use the default
     // spatial discretization scheme, which is not accurate for certain
@@ -211,8 +210,8 @@ int main() {
     // For coordinate independent velocity functions
     // this numerical scheme is superior though.
     // However, it is slower.
-    // advectionKernel.setSpatialScheme(
-    //     ls::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
+    advectionKernel.setSpatialScheme(
+        ls::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
   }
 
   // advect the level set until 50s have passed

diff --git a/include/viennals/lsAdvect.hpp b/include/viennals/lsAdvect.hpp
@@ -24,7 +24,7 @@
 #include <lsLocalLaxFriedrichsAnalytical.hpp>
 #include <lsLocalLocalLaxFriedrichs.hpp>
 #include <lsStencilLocalLaxFriedrichsScalar.hpp>
-#include <lsWENO5.hpp>
+#include <lsWENO.hpp>
 
 // Include implementation of time integration schemes
 #include <lsAdvectIntegrationSchemes.hpp>
@@ -147,12 +147,12 @@ template <class T, int D> class Advect {
             // move neighborIterator to current position
             neighborIterator.goToIndicesSequential(indices);
 
-            Vec3D<T> coords;
+            Vec3D<T> coords{};
             for (unsigned i = 0; i < D; ++i) {
               coords[i] = indices[i] * gridDelta;
             }
 
-            Vec3D<T> normal = {};
+            Vec3D<T> normal{};
             T normalModulus = 0.;
             for (unsigned i = 0; i < D; ++i) {
               const T phiPos = neighborIterator.getNeighbor(i).getValue();
@@ -664,10 +664,15 @@ template <class T, int D> class Advect {
       auto is = lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
+      // Instantiate WENO with order 3
+      auto is = lsInternal::WENO<T, D, 3>(levelSets.back(), velocities,
+                                          dissipationAlpha);
+      currentTimeStep = integrateTime(is, maxTimeStep);
     } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
-      // Instantiate WENO5 with order 3 (neighbors +/- 3)
-      auto is = lsInternal::WENO5<T, D, 3>(levelSets.back(), velocities,
-                                           dissipationAlpha);
+      // Instantiate WENO with order 5
+      auto is = lsInternal::WENO<T, D, 5>(levelSets.back(), velocities,
+                                          dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
     } else {
       VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
@@ -998,9 +1003,10 @@ template <class T, int D> class Advect {
                SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>::prepareLS(
           levelSets.back());
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
+      lsInternal::WENO<T, D, 3>::prepareLS(levelSets.back());
     } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
-      // WENO5 requires a stencil radius of 3 (template parameter 3)
-      lsInternal::WENO5<T, D, 3>::prepareLS(levelSets.back());
+      lsInternal::WENO<T, D, 5>::prepareLS(levelSets.back());
     } else {
       VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
     }

diff --git a/include/viennals/lsAdvectIntegrationSchemes.hpp b/include/viennals/lsAdvectIntegrationSchemes.hpp
@@ -18,7 +18,8 @@ enum class SpatialSchemeEnum : unsigned {
   LOCAL_LAX_FRIEDRICHS_1ST_ORDER = 7,
   LOCAL_LAX_FRIEDRICHS_2ND_ORDER = 8,
   STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER = 9,
-  WENO_5TH_ORDER = 10
+  WENO_3RD_ORDER = 10,
+  WENO_5TH_ORDER = 11
 };
 
 // Legacy naming (deprecated, will be removed in future versions)

diff --git a/include/viennals/lsCalculateNormalVectors.hpp b/include/viennals/lsCalculateNormalVectors.hpp
@@ -134,12 +134,12 @@ template <class T, int D> class CalculateNormalVectors {
           continue;
         } else if (std::abs(center.getValue()) > maxValue) {
           // push an empty vector to keep ordering correct
-          Vec3D<T> tmp = {};
+          Vec3D<T> tmp{};
           normalVectors.push_back(tmp);
           continue;
         }
 
-        Vec3D<T> n;
+        Vec3D<T> n{};
 
         T denominator = 0;
         for (int i = 0; i < D; i++) {

diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp
@@ -41,8 +41,8 @@ template <class T, int D> class CalculateVisibilities {
       std::fill(visibilities.begin(), visibilities.end(), static_cast<T>(1.0));
     } else {
       // *** Determine extents of domain ***
-      Vec3D<T> minDefinedPoint;
-      Vec3D<T> maxDefinedPoint;
+      Vec3D<T> minDefinedPoint{};
+      Vec3D<T> maxDefinedPoint{};
       // Initialize with extreme values
       for (int i = 0; i < D; ++i) {
         minDefinedPoint[i] = std::numeric_limits<T>::max();
@@ -87,7 +87,7 @@ template <class T, int D> class CalculateVisibilities {
             continue;
 
           // Starting position of the point
-          Vec3D<T> currentPos;
+          Vec3D<T> currentPos{};
           for (int i = 0; i < D; ++i) {
             currentPos[i] = it.getStartIndices(i);
           }

diff --git a/include/viennals/lsCompareCriticalDimensions.hpp b/include/viennals/lsCompareCriticalDimensions.hpp
@@ -339,7 +339,8 @@ template <class T, int D = 2> class CompareCriticalDimensions {
         continue;
 
       // Create points for target and sample positions
-      Vec3D<T> coordTarget = {0.0, 0.0, 0.0}, coordSample = {0.0, 0.0, 0.0};
+      Vec3D<T> coordTarget{};
+      Vec3D<T> coordSample{};
 
       for (int i = 0; i < D; ++i) {
         if (i == result.measureDimension) {

diff --git a/include/viennals/lsDomain.hpp b/include/viennals/lsDomain.hpp
@@ -2,6 +2,7 @@
 
 #include <lsPreCompileMacros.hpp>
 
+#include <cstring>
 #include <limits>
 
 #include <hrleDomain.hpp>
@@ -221,7 +222,7 @@ template <class T, int D> class Domain {
     // Check identifier
     char identifier[8];
     stream.read(identifier, 8);
-    if (std::string(identifier).compare(0, 8, "lsDomain")) {
+    if (std::memcmp(identifier, "lsDomain", 8) != 0) {
       Logger::getInstance()
           .addError(
               "Reading Domain from stream failed. Header could not be found.")

diff --git a/include/viennals/lsEngquistOsher.hpp b/include/viennals/lsEngquistOsher.hpp
@@ -121,7 +121,7 @@ template <class T, int D, int order> class EngquistOsher {
 
     // Calculate normal vector for velocity calculation
     // use std::array since it will be exposed to interface
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     if (calculateNormalVectors) {
       T denominator = 0;
       for (int i = 0; i < D; i++) {

diff --git a/include/viennals/lsExtrude.hpp b/include/viennals/lsExtrude.hpp
@@ -19,7 +19,7 @@ template <class T> class Extrude {
   SmartPointer<Domain<T, 3>> outputLevelSet = nullptr;
   Vec2D<T> extent = {0., 0.};
   int extrusionAxis = 0;
-  std::array<BoundaryConditionEnum, 3> boundaryConds = {};
+  std::array<BoundaryConditionEnum, 3> boundaryConds{};
 
 public:
   Extrude() = default;

diff --git a/include/viennals/lsGeometricAdvect.hpp b/include/viennals/lsGeometricAdvect.hpp
@@ -482,7 +482,7 @@ template <class T, int D> class GeometricAdvect {
     {
       std::vector<T> scalarData;
       for (auto it = newPoints[0].begin(); it != newPoints[0].end(); ++it) {
-        Vec3D<T> node{0., 0., 0.};
+        Vec3D<T> node{};
         for (unsigned i = 0; i < D; ++i) {
           node[i] = T((it->first)[i]) * gridDelta;
         }

diff --git a/include/viennals/lsGeometricAdvectDistributions.hpp b/include/viennals/lsGeometricAdvectDistributions.hpp
@@ -108,7 +108,7 @@ class SphereDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -radius;
       bounds[2 * i + 1] = radius;
@@ -157,7 +157,7 @@ class BoxDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -posExtent[i];
       bounds[2 * i + 1] = posExtent[i];
@@ -240,7 +240,7 @@ class CustomSphereDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -maxRadius_;
       bounds[2 * i + 1] = maxRadius_;

diff --git a/include/viennals/lsLaxFriedrichs.hpp b/include/viennals/lsLaxFriedrichs.hpp
@@ -61,7 +61,7 @@ template <class T, int D, int order> class LaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions

diff --git a/include/viennals/lsLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLaxFriedrichs.hpp
@@ -86,7 +86,7 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
@@ -183,7 +183,7 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
     }
 
     // calculate alphas
-    T alpha[D] = {};
+    T alpha[D]{};
     {
       // alpha calculation is always on order 1 stencil
       const viennahrle::IndexType minIndex = -1;
@@ -192,11 +192,11 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
 
       viennahrle::Index<D> neighborIndex(minIndex);
       for (unsigned i = 0; i < numNeighbors; ++i) {
-        Vec3D<T> coords;
+        Vec3D<T> coords{};
         for (unsigned dir = 0; dir < D; ++dir) {
           coords[dir] = coordinate[dir] + neighborIndex[dir] * gridDelta;
         }
-        Vec3D<T> normal = {};
+        Vec3D<T> normal{};
         double normalModulus = 0.;
         auto center = neighborIterator.getNeighbor(neighborIndex).getValue();
         for (unsigned dir = 0; dir < D; ++dir) {

diff --git a/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp b/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
@@ -85,7 +85,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
@@ -182,7 +182,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
     }
 
     // calculate alphas
-    T alpha[D] = {};
+    T alpha[D]{};
     {
       // alpha calculation is always on order 1 stencil
       constexpr viennahrle::IndexType minIndex = -1;
@@ -192,7 +192,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
       viennahrle::Index<D> neighborIndex(minIndex);
       for (unsigned i = 0; i < numNeighbors; ++i) {
 
-        Vec3D<T> normal = {};
+        Vec3D<T> normal{};
         auto center = neighborIterator.getNeighbor(neighborIndex).getValue();
         for (unsigned dir = 0; dir < D; ++dir) {
           viennahrle::Index<D> unity(0);

diff --git a/include/viennals/lsLocalLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLocalLaxFriedrichs.hpp
@@ -64,7 +64,7 @@ template <class T, int D, int order> class LocalLocalLaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions

diff --git a/include/viennals/lsMakeGeometry.hpp b/include/viennals/lsMakeGeometry.hpp
@@ -441,10 +441,9 @@ template <class T, int D> class MakeGeometry {
 
         // create and insert points at base
         for (double angle = 0.; angle < limit; angle += smallAngle) {
-          Vec3D<T> point;
+          Vec3D<T> point{};
           point[0] = cylinder->radius * std::cos(angle);
           point[1] = cylinder->radius * std::sin(angle);
-          point[2] = 0.0;
           points.push_back(point);
           mesh->insertNextNode(point);
         }

diff --git a/include/viennals/lsMesh.hpp b/include/viennals/lsMesh.hpp
@@ -29,8 +29,8 @@ template <class T = double> class Mesh {
   std::vector<std::array<unsigned, 8>> hexas;
   PointData<T> pointData;
   PointData<T> cellData;
-  Vec3D<T> minimumExtent;
-  Vec3D<T> maximumExtent;
+  Vec3D<T> minimumExtent{};
+  Vec3D<T> maximumExtent{};
 
 private:
   // iterator typedef
@@ -272,6 +272,8 @@ template <class T = double> class Mesh {
     hexas.clear();
     pointData.clear();
     cellData.clear();
+    minimumExtent = Vec3D<T>{};
+    maximumExtent = Vec3D<T>{};
   }
 
   void print() {

diff --git a/include/viennals/lsPointData.hpp b/include/viennals/lsPointData.hpp
@@ -3,6 +3,7 @@
 #include <lsPreCompileMacros.hpp>
 
 #include <cstdint>
+#include <cstring>
 #include <utility>
 #include <vector>
 
@@ -385,7 +386,7 @@ class PointData {
   std::istream &deserialize(std::istream &stream) {
     char identifier[11];
     stream.read(identifier, 11);
-    if (std::string(identifier).compare(0, 11, "lsPointData")) {
+    if (std::memcmp(identifier, "lsPointData", 11) != 0) {
       Logger::getInstance()
           .addError("Reading PointData from stream failed. Header could "
                     "not be found.")