diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3ea95a65f7e0d937949b8c5a359cbc47d712ba0c..1ce98ccfbbcdb556a308777ece789c40caa03635 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -54,7 +54,7 @@ if(ENABLE_GPU)
elseif ( CUDA_VERSION_MAJOR EQUAL 10 AND CUDA_VERSION_MINOR EQUAL 2 )
message("CUDA is compatible")
set(WARNING_SUPPRESSION_AND_OPTION_NVCC -Xcudafe "--display_error_number --diag_suppress=2976 --diag_suppress=2977 --diag_suppress=2979 --diag_suppress=186" --expt-extended-lambda)
-elseif ( CUDA_VERSION_MAJOR EQUAL 11 AND CUDA_VERSION_MINOR EQUAL 0 )
+ elseif ( CUDA_VERSION_MAJOR EQUAL 11 AND CUDA_VERSION_MINOR EQUAL 0 )
message("CUDA is compatible")
set(WARNING_SUPPRESSION_AND_OPTION_NVCC -Xcudafe "--display_error_number --diag_suppress=3059 --diag_suppress=3058 --diag_suppress=3057 --diag_suppress=3056 --diag_suppress=611 --diag_suppress=186" --expt-extended-lambda)
else()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4371f488abd5b16454d4357e3fa5dccac4022115..a20a228da68da0cb7f5d8a4e0ecfe6d8c629ca27 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -143,6 +143,8 @@ install(FILES Operators/Vector/vector_dist_operators_extensions.hpp
install(FILES Draw/DrawParticles.hpp
Draw/PointIterator.hpp
Draw/PointIteratorSkin.hpp
+ Draw/DrawCircle.hpp
+ Draw/DrawSphere.hpp
DESTINATION openfpm_numerics/include/Draw )
install(FILES interpolation/interpolation.hpp
@@ -150,6 +152,14 @@ install(FILES interpolation/interpolation.hpp
interpolation/z_spline.hpp
DESTINATION openfpm_numerics/include/interpolation )
+install(FILES level_set/redistancing_Sussman/ComputeGradient.hpp
+ level_set/redistancing_Sussman/GaussFilter.hpp
+ level_set/redistancing_Sussman/HelpFunctions.hpp
+ level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
+ level_set/redistancing_Sussman/NarrowBand.hpp
+ level_set/redistancing_Sussman/RedistancingSussman.hpp
+ DESTINATION openfpm_numerics/include/level_set/redistancing_Sussman)
+
install(FILES DMatrix/EMatrix.hpp
DESTINATION openfpm_numerics/include/DMatrix )
diff --git a/src/Draw/DrawCircle.hpp b/src/Draw/DrawCircle.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..062e4ca5f6167ce51c1cc8d810cb592183884242
--- /dev/null
+++ b/src/Draw/DrawCircle.hpp
@@ -0,0 +1,80 @@
+//
+// Created by jstark on 2020-08-18.
+//
+/**
+ * @file Circle.hpp
+ *
+ * @brief Header file containing functions for creating a filled 2D circle of defined radius.
+ *
+ * @author Justina Stark
+ * @date August 2020
+ */
+#ifndef ACCURACY_TESTS_CIRCLE_HPP
+#define ACCURACY_TESTS_CIRCLE_HPP
+
+#include <iostream>
+
+/**@brief Checks if a point lays inside a circle of given radius and given center coordinates.
+ *
+ * @tparam coord_type Inferred type of the point coordinates x and y.
+ * @tparam radius_type Inferred type of radius the circle is supposed to have.
+ * @param x_coord X-coordinate of the point for that laying inside the circle should be checked.
+ * @param y_coord Y-coordinate of the point for that laying inside the circle should be checked.
+ * @param radius Radius of the filled circle.
+ * @param center_x X-coordinate of the circle center.
+ * @param center_y Y-coordinate of the circle center.
+ * @return True, if point lays inside the circle; false, if outside.
+ */
+template <typename coord_type, typename radius_type>
+bool inside_circle(coord_type x_coord, coord_type y_coord, radius_type radius, float center_x=0, float center_y=0)
+{
+ return (x_coord - center_x) * (x_coord - center_x) + (y_coord - center_y) * (y_coord - center_y) <= radius * radius;
+}
+
+/**@brief Creates a filled circle of given radius an given center coordinates on a 2D OpenFPM grid.
+ *
+ * @details The circle is represented in one of the grid properties via the following indicator function:
+ * @f[ \phi_{\text{indicator}} = \begin{cases}
+ * +1 & \text{point lies inside the circle} \\
+ * -1 & \text{point lies outside the circle} \\
+ * \end{cases} @f]
+ *
+ * @tparam Phi_0 Index of grid property that should store the indicator function representing the filled circle.
+ * @tparam grid_type Inferred type of the input grid.
+ * @tparam radius_type Inferred type of radius the circle is supposed to have.
+ * @param grid Input OpenFPM grid.
+ * @param radius Radius of the filled circle.
+ * @param center_x X-coordinate of the circle center.
+ * @param center_y Y-coordinate of the circle center.
+ */
+template <size_t Phi_0, typename grid_type, typename radius_type>
+void init_grid_with_circle(grid_type & grid, radius_type radius, float center_x=0, float center_y=0)
+{
+ const size_t x = 0;
+ const size_t y = 1;
+ // assign pixel values to domain. For each pixel get factor_refinement number of grid points with corresponding value
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto gkey = grid.getGKey(key);
+ auto spacing = grid.getSpacing();
+ auto i = gkey.get(x);
+ auto j = gkey.get(y);
+ double x_coord = i * spacing[x];
+ double y_coord = j * spacing[y];
+
+ if (inside_circle(x_coord, y_coord, radius, center_x, center_y))
+ {
+ grid.template get<Phi_0> (key) = 1;
+ }
+ else
+ {
+ grid.template get<Phi_0> (key) = -1;
+ }
+ ++dom;
+ }
+}
+
+
+#endif //ACCURACY_TESTS_CIRCLE_HPP
diff --git a/src/Draw/DrawSphere.hpp b/src/Draw/DrawSphere.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..576bd17a78d06a0e1e73fd15f8e9157150efa7a7
--- /dev/null
+++ b/src/Draw/DrawSphere.hpp
@@ -0,0 +1,85 @@
+//
+// Created by jstark on 2020-05-17.
+//
+/**
+ * @file Sphere.hpp
+ *
+ * @brief Header file containing functions for creating a filled 3D sphere of defined radius.
+ *
+ * @author Justina Stark
+ * @date May 2020
+ */
+#ifndef ACCURACY_TESTS_SPHERE_HPP
+#define ACCURACY_TESTS_SPHERE_HPP
+
+#include <iostream>
+/**@brief Checks if a point lays inside a sphere of given radius and given center coordinates.
+ *
+ * @tparam coord_type Inferred type of the point coordinates x, y and z.
+ * @tparam radius_type Inferred type of radius the sphere is supposed to have.
+ * @param x_coord X-coordinate of the point for that laying inside the sphere should be checked.
+ * @param y_coord Y-coordinate of the point for that laying inside the sphere should be checked.
+ * @param z_coord Z-coordinate of the point for that laying inside the sphere should be checked.
+ * @param radius Radius of the filled sphere.
+ * @param center_x X-coordinate of the sphere center.
+ * @param center_y Y-coordinate of the sphere center.
+ * @param center_z Z-coordinate of the sphere center.
+ * @return True, if point lays inside the sphere; false, if outside.
+ */
+template <typename coord_type, typename radius_type>
+bool inside_sphere(coord_type x_coord, coord_type y_coord, coord_type z_coord, radius_type radius, float center_x=0, float center_y=0, float center_z=0)
+{
+ return (x_coord - center_x) * (x_coord - center_x) + (y_coord - center_y) * (y_coord - center_y) + (z_coord - center_z) * (z_coord - center_z) <= radius * radius;
+}
+
+/**@brief Creates a filled sphere of given radius an given center coordinates on a 3D OpenFPM grid.
+ *
+ * @details The sphere is represented in one of the grid properties via the following indicator function:
+ * @f[ \phi_{\text{indicator}} = \begin{cases}
+ * +1 & \text{point lies inside the sphere} \\
+ * -1 & \text{point lies outside the sphere} \\
+ * \end{cases} @f]
+ *
+ * @tparam Phi_0 Index of grid property that should store the indicator function representing the filled sphere.
+ * @tparam grid_type Inferred type of the input grid.
+ * @tparam radius_type Inferred type of radius the sphere is supposed to have.
+ * @param grid Input OpenFPM grid.
+ * @param radius Radius of the filled sphere.
+ * @param center_x X-coordinate of the sphere center.
+ * @param center_y Y-coordinate of the sphere center.
+ * @param center_z Z-coordinate of the sphere center.
+ */
+template <size_t Phi_0, typename grid_type, typename radius_type>
+void init_grid_with_sphere(grid_type & grid, radius_type radius, float center_x=0, float center_y=0, float center_z=0)
+{
+ const size_t x = 0;
+ const size_t y = 1;
+ const size_t z = 2;
+ // assign pixel values to domain. For each pixel get factor_refinement number of grid points with corresponding value
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto gkey = grid.getGKey(key);
+ auto spacing = grid.getSpacing();
+ auto i = gkey.get(x);
+ auto j = gkey.get(y);
+ auto k = gkey.get(z);
+ double x_coord = i * spacing[x];
+ double y_coord = j * spacing[y];
+ double z_coord = k * spacing[z];
+
+ if (inside_sphere(x_coord, y_coord, z_coord, radius, center_x, center_y, center_z))
+ {
+ grid.template get<Phi_0> (key) = 1;
+ }
+ else
+ {
+ grid.template get<Phi_0> (key) = -1;
+ }
+ ++dom;
+ }
+}
+
+
+#endif //ACCURACY_TESTS_SPHERE_HPP
diff --git a/src/FiniteDifference/Eno_Weno.hpp b/src/FiniteDifference/Eno_Weno.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc2bf99a0d40e597643f09971e0bfe8364760adb
--- /dev/null
+++ b/src/FiniteDifference/Eno_Weno.hpp
@@ -0,0 +1,153 @@
+/*
+ *
+ * Based on the description in Osher and Fedkiw, Level Set Methods and Dynamic Implicit Surfaces
+ * For HJ ENO method - Sec. 3.3
+ * For HJ WENO method - Sec. 3.4
+*/
+
+#ifndef __ENO_WENO_HPP__
+#define __ENO_WENO_HPP__
+
+#include "Grid/grid_dist_id.hpp"
+
+double adjustWeights(double v1, double v2, double v3, double v4, double v5)
+{
+ double phix1, phix2, phix3;
+ double s1, s2, s3;
+ double a1, a2, a3;
+ double eps = 1e-6;
+ double w1, w2, w3;
+
+ // Eqs. (3.25), (3.26), (3.27)
+ phix1 = (v1 / 3.0) - (7.0 * v2 / 6.0) + (11.0 * v3 / 6.0);
+ phix2 = -(v2 / 6.0) + (5.0 * v3 / 6.0) + (v4 / 3.0);
+ phix3 = (v3 / 3.0) + (5.0 * v4 / 6.0) - (v5 / 6.0);
+
+ // Eqs. (3.32), (3.33), (3.34)
+ s1 = (13.0 / 12.0) * (v1 - 2*v2 + v3) * (v1 - 2*v2 + v3)
+ + 0.25 * (v1 - 4*v2 + 3*v3) * (v1 - 4*v2 + 3*v3);
+ s2 = (13.0 / 12.0) * (v2 - 2*v3 + v4) * (v2 - 2*v3 + v4)
+ + 0.25 * (v2 - v4) * (v2 - v4);
+ s3 = (13.0 / 12.0) * (v3 - 2*v4 + v5) * (v3 - 2*v4 + v5)
+ + 0.25 * (3*v3 - 4*v4 + v5) * (3*v3 - 4*v4 + v5);
+
+ // Eqs. (3.35), (3.36), (3.37)
+ a1 = 0.1 / ((s1 + eps)*(s1 + eps));
+ a2 = 0.6 / ((s2 + eps)*(s2 + eps));
+ a3 = 0.3 / ((s3 + eps)*(s3 + eps));
+
+ // Eqs. (3.39), (3.40), (3.41)
+ w1 = a1 / (a1 + a2 + a3);
+ w2 = a2 / (a1 + a2 + a3);
+ w3 = a3 / (a1 + a2 + a3);
+
+ // Eq. (3.28)
+ return (w1 * phix1 + w2 * phix2 + w3 * phix3);
+}
+
+template<typename GridDist, typename GridKey, unsigned int field, unsigned int x>
+double WENO_5_Plus(GridKey key, GridDist &gd)
+{
+ double c2, c3, c4, c5, c6;
+ double coeff = 1.0 / gd.spacing(x);
+ double df;
+
+ c2 = (gd.template get<field>(key.move(x, -1)) - gd.template get<field>(key.move(x, -2))) * coeff;
+ c3 = (gd.template get<field>(key) - gd.template get<field>(key.move(x, -1))) * coeff;
+ c4 = (gd.template get<field>(key.move(x,1)) - gd.template get<field>(key)) * coeff;
+ c5 = (gd.template get<field>(key.move(x, 2)) - gd.template get<field>(key.move(x, 1))) * coeff;
+ c6 = (gd.template get<field>(key.move(x, 3)) - gd.template get<field>(key.move(x, 2))) * coeff;
+ df = adjustWeights(c6, c5, c4, c3, c2);
+
+ return df;
+}
+
+template<typename GridDist, typename GridKey, unsigned int field, unsigned int x>
+double WENO_5_Minus(GridKey key, GridDist &gd)
+{
+ double c1, c2, c3, c4, c5;
+ double coeff = 1.0 / gd.spacing(x);
+ double df;
+
+ c1 = (gd.template get<field>(key.move(x, -2)) - gd.template get<field>(key.move(x, -3))) * coeff;
+ c2 = (gd.template get<field>(key.move(x, -1)) - gd.template get<field>(key.move(x, -2))) * coeff;
+ c3 = (gd.template get<field>(key) - gd.template get<field>(key.move(x, -1))) * coeff;
+ c4 = (gd.template get<field>(key.move(x,1)) - gd.template get<field>(key)) * coeff;
+ c5 = (gd.template get<field>(key.move(x, 2)) - gd.template get<field>(key.move(x, 1))) * coeff;
+
+ df = adjustWeights(c1, c2, c3, c4, c5);
+
+ return df;
+}
+
+template<typename GridDist, typename GridKey, unsigned int field, unsigned int x>
+double ENO_3_Plus(GridKey key, GridDist &gd)
+{
+ double q1x, q2x, q3x;
+ double gridsize = gd.spacing(x);
+ double coeff = 1.0 / gridsize;
+
+ q1x = (gd.template get<field>(key.move(x,1)) - gd.template get<field>(key)) * coeff;
+ double d2ix = (gd.template get<field>(key.move(x,1)) - 2*gd.template get<field>(key) + gd.template get<field>(key.move(x,-1))) * 0.5 * coeff * coeff;
+ double d2ip1x = (gd.template get<field>(key.move(x,2)) - 2*gd.template get<field>(key.move(x,1)) + gd.template get<field>(key)) * 0.5 * coeff * coeff;
+ if(abs(d2ix) <= abs(d2ip1x))
+ {
+ q2x = -d2ix * gridsize;
+ double d3p1hx = (gd.template get<field>(key.move(x,1)) - 3*gd.template get<field>(key) + 3*gd.template get<field>(key.move(x,-1)) - gd.template get<field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
+ double d3p3hx = (gd.template get<field>(key.move(x,2)) - 3*gd.template get<field>(key.move(x,1)) + 3*gd.template get<field>(key) - gd.template get<field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
+ if(abs(d3p1hx) <= abs(d3p3hx))
+ q3x = -d3p1hx * gridsize * gridsize;
+ else
+ q3x = -d3p3hx * gridsize * gridsize;
+ }
+ else
+ {
+ q2x = -d2ip1x * gridsize;
+ double d3p1hx = (gd.template get<field>(key.move(x,2)) - 3*gd.template get<field>(key.move(x,1)) + 3*gd.template get<field>(key) - gd.template get<field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
+ double d3p3hx = (gd.template get<field>(key.move(x,3)) - 3*gd.template get<field>(key.move(x,2)) + 3*gd.template get<field>(key.move(x,1)) - gd.template get<field>(key)) * coeff * coeff * coeff / 6.0;
+ if(abs(d3p1hx) <= abs(d3p3hx))
+ q3x = d3p1hx * 2 * gridsize * gridsize;
+ else
+ q3x = d3p3hx * 2 * gridsize * gridsize;
+
+ }
+ return (q1x + q2x + q3x);
+}
+
+template<typename GridDist, typename GridKey, unsigned int field, unsigned int x>
+double ENO_3_Minus(GridKey key, GridDist &gd)
+{
+ double q1x, q2x, q3x;
+ double gridsize = gd.spacing(x);
+ double coeff = 1.0 / gridsize;
+
+ q1x = (gd.template get<field>(key) - gd.template get<field>(key.move(x, -1))) * coeff;
+ double d2im1x = (gd.template get<field>(key) - 2*gd.template get<field>(key.move(x, -1)) + gd.template get<field>(key.move(x, -2))) * 0.5 * coeff * coeff;
+ double d2ix = (gd.template get<field>(key.move(x,1)) - 2*gd.template get<field>(key) + gd.template get<field>(key.move(x,-1))) * 0.5 * coeff * coeff;
+
+ if(abs(d2im1x) <= abs(d2ix))
+ {
+ q2x = d2im1x * gridsize;
+ double d3p1hx = (gd.template get<field>(key) - 3*gd.template get<field>(key.move(x,-1)) + 3*gd.template get<field>(key.move(x,-2)) - gd.template get<field>(key.move(x,-3))) * coeff * coeff * coeff / 6.0;
+ double d3p3hx = (gd.template get<field>(key.move(x,1)) - 3*gd.template get<field>(key) + 3*gd.template get<field>(key.move(x,-1)) - gd.template get<field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
+ if(abs(d3p1hx) <= abs(d3p3hx))
+ q3x = d3p1hx * 2 * gridsize * gridsize;
+ else
+ q3x = d3p3hx * 2 * gridsize * gridsize;
+ }
+ else
+ {
+ q2x = d2ix * gridsize;
+ double d3p1hx = (gd.template get<field>(key.move(x,1)) - 3*gd.template get<field>(key) + 3*gd.template get<field>(key.move(x,-1)) - gd.template get<field>(key.move(x,-2))) * coeff * coeff * coeff / 6.0;
+ double d3p3hx = (gd.template get<field>(key.move(x,2)) - 3*gd.template get<field>(key.move(x,1)) + 3*gd.template get<field>(key) - gd.template get<field>(key.move(x,-1))) * coeff * coeff * coeff / 6.0;
+ if(abs(d3p1hx) <= abs(d3p3hx))
+ q3x = -d3p1hx * gridsize * gridsize;
+ else
+ q3x = -d3p3hx * gridsize * gridsize;
+ }
+
+ return (q1x + q2x + q3x);
+}
+
+#endif
+
diff --git a/src/level_set/redistancing_Sussman/ComputeGradient.hpp b/src/level_set/redistancing_Sussman/ComputeGradient.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..21ef4431c9c1d1af724732848dfbca35d4245f5b
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/ComputeGradient.hpp
@@ -0,0 +1,184 @@
+//
+// Created by jstark on 2020-05-08.
+//
+/**
+ * @file ComputeGradient.hpp
+ *
+ * @brief Header file containing the functions to compute an upwind gradient and the magnitude of gradient.
+ *
+ * @details Upwind gradient is computed according to <a href="https://www.researchgate.net/publication/2642502_An_Efficient_Interface_Preserving_Level_Set_Re
+ * -Distancing_Algorithm_And_Its_Application_To_Interfacial_Incompressible_Fluid_Flow.html">M. Sussman and E. Fatemi,
+ * “Efficient, interface-preserving level set redistancing algorithm and its application to interfacial
+ * incompressible fluid flow†(1999)</a> ).
+ *
+ * @author Justina Stark
+ * @date May 2020
+ */
+#ifndef REDISTANCING_SUSSMAN_COMPUTEGRADIENT_HPP
+#define REDISTANCING_SUSSMAN_COMPUTEGRADIENT_HPP
+
+// Include standard library header files
+#include <iostream>
+#include <fstream>
+#include <typeinfo>
+#include <cmath>
+#include <array>
+// Include OpenFPM header files
+#include "Grid/grid_dist_id.hpp"
+
+/**@brief Computes the forward finite difference of a scalar property on the current grid node.
+ *
+ * @tparam Phi Size_t index of property for which the gradient should be computed.
+ * @tparam gridtype Type of input grid.
+ * @tparam keytype Type of key variable.
+ * @param grid Grid, on which the gradient should be computed.
+ * @param key Key that contains the index of the current grid node.
+ * @param d Variable (size_t) that contains the dimension.
+ * @return Forward finite difference for the property under index Phi on the current node with index key.
+ */
+template <size_t Phi, typename gridtype, typename keytype>
+double forward_FD(gridtype & grid, keytype & key, size_t d)
+{
+ return (grid.template get<Phi> (key.move(d, 1)) - grid.template get<Phi> (key)) / grid.getSpacing()[d];
+}
+/**@brief Computes the backward finite difference of a scalar property on the current grid node.
+ *
+ * @tparam Phi Size_t index of property for which the gradient should be computed.
+ * @tparam gridtype Type of input grid.
+ * @tparam keytype Type of key variable.
+ * @param grid Grid, on which the gradient should be computed.
+ * @param key Key that contains the index of the current grid node.
+ * @param d Variable (size_t) that contains the dimension.
+ * @return Backward finite difference for the property under index Phi on the current node with index key.
+ */
+template <size_t Phi, typename gridtype, typename keytype>
+double backward_FD(gridtype & grid, keytype & key, size_t d)
+{
+ return (grid.template get<Phi> (key) - grid.template get<Phi> (key.move(d, -1))) / grid.getSpacing()[d];
+}
+/**@brief Get the first order upwind finite difference of a scalar property on the current grid node.
+ *
+ * @details Calls #forward_FD() and #backward_FD() and chooses between dplus and dminus by upwinding.
+ *
+ * @tparam Phi Size_t index of property for which the gradient should be computed.
+ * @tparam Phi_0_sign Size_t index of property that contains the initial sign of that property for which the upwind FD
+ * should be computed.
+ * @tparam gridtype Type of input grid.
+ * @tparam keytype Type of key variable.
+ * @param grid Grid, on which the gradient should be computed.
+ * @param key Key that contains the index of the current grid node.
+ * @param d Variable (size_t) that contains the dimension.
+ * @return First order upwind finite difference for the property under index Phi on the current node with index key.
+ */
+template <size_t Phi, size_t Phi_0_sign, typename gridtype, typename keytype>
+double upwind_FD(gridtype &grid, keytype &key, size_t d)
+{
+ double dplus = 0, dminus = 0, dphi = 0;
+
+ dplus = forward_FD<Phi>(grid, key, d);
+ dminus = backward_FD<Phi>(grid, key, d);
+ int sign_phi0 = grid.template get<Phi_0_sign> (key);
+
+ if (dplus * sign_phi0 < 0
+ && (dminus + dplus) * sign_phi0 < 0) dphi = dplus;
+
+ else if (dminus * sign_phi0 > 0
+ && (dminus + dplus) * sign_phi0 > 0) dphi = dminus;
+
+ else if (dminus * sign_phi0 < 0
+ && dplus * sign_phi0 > 0) dphi = 0;
+
+ //else if (-10e-12 <= dminus <= 10e-12 && -10e-12 <= dplus <= 10e-12) dphi = 0;
+
+ return dphi;
+}
+/**@brief Computes 1st order upwind finite difference inside and forward/backward finite difference at the grid borders.
+ *
+ * @details Checks if a point lays within the grid or at the border. For the internal grid points it calls upwind
+ * finite difference #upwind_FD(). For the border points, simply #forward_FD() and #backward_FD() is used,
+ * respectively,
+ * depending on the side of the border.
+ *
+ * @tparam Phi Size_t index of property for which the gradient should be computed.
+ * @tparam Phi_0_sign Size_t index of property that contains the initial sign of that property for which the upwind FD
+ * should be computed.
+ * @tparam Phi_grad Size_t index of property where the gradient result should be stored.
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the gradient should be computed.
+ */
+// Use upwinding for inner grid points and one sided backward / forward stencil at border grid points
+template <size_t Phi, size_t Phi_0_sign, size_t Phi_grad, typename gridtype>
+void get_first_order_gradient_depending_on_point_position(gridtype &grid)
+{
+ grid.template ghost_get<Phi>();
+ auto dom = grid.getDomainIterator();
+ while (dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_g = grid.getGKey(key);
+
+ for(size_t d = 0; d < gridtype::dims; d++ )
+ {
+ if (key_g.get(d) > 0 && key_g.get(d) < grid.size(d) - 1) // if point lays not at the border of the grid
+ {
+ grid.template get<Phi_grad> (key) [d] = upwind_FD<Phi, Phi_0_sign>(grid, key, d);
+ }
+
+ else if (key_g.get(d) == 0) // if point lays at left border, use right sided kernel
+ {
+ grid.template get<Phi_grad> (key) [d] = forward_FD<Phi>(grid, key, d);
+ }
+
+ else if (key_g.get(d) >= grid.size(d) - 1) // if point lays at right border, use left sided kernel
+ {
+ grid.template get<Phi_grad> (key) [d] = backward_FD<Phi>(grid, key, d);
+ }
+ }
+ ++dom;
+ }
+}
+/**@brief Calls #get_first_order_gradient_depending_on_point_position(). Computes 1st order upwind finite difference.
+ *
+* @tparam Phi Size_t index of property for which the gradient should be computed.
+ * @tparam Phi_0_sign Size_t index of property that contains the initial sign of that property for which the upwind FD
+ * should be computed.
+ * @tparam Phi_grad Size_t index of property where the gradient result should be stored.
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the gradient should be computed.
+ */
+/*Approximate the initial gradients between neighbor particles using fwd and bwd differences.
+ These gradients are needed for the re-distancing step according to Sussman & Fatemi 1999.*/
+template <size_t Phi_in, size_t Phi_0_sign, size_t Phi_grad_out, typename gridtype>
+void get_upwind_gradient(gridtype & grid)
+{
+ get_first_order_gradient_depending_on_point_position<Phi_in, Phi_0_sign, Phi_grad_out>(grid);
+}
+
+/**@brief Computes the magnitude of the gradient (L2-norm of gradient vector).
+ *
+ * @tparam Phi_grad_in Size_t index of property that contains the gradient.
+ * @tparam Phi_magnOfGrad_out Size_t index of property where the magnitude of gradient should be stored.
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the magnitude of gradient should be computed.
+ */
+template <size_t Phi_grad_in, size_t Phi_magnOfGrad_out, typename gridtype>
+void get_gradient_magnitude(gridtype & grid)
+{
+ grid.template ghost_get<Phi_grad_in>();
+ auto dom = grid.getDomainGhostIterator();
+ while(dom.isNext())
+ {
+ double sum = 0;
+ auto key = dom.get();
+ for(size_t d = 0; d < gridtype::dims; d++)
+ {
+ sum += grid.template get<Phi_grad_in> (key)[d] * grid.template get<Phi_grad_in> (key)[d];
+ }
+ grid.template get<Phi_magnOfGrad_out> (key) = sqrt(sum);
+ ++dom;
+ }
+}
+
+
+
+#endif //REDISTANCING_SUSSMAN_COMPUTEGRADIENT_HPP
diff --git a/src/level_set/redistancing_Sussman/GaussFilter.hpp b/src/level_set/redistancing_Sussman/GaussFilter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4121f634ed924aaf32c2731081c9cc6ba2f3dc4d
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/GaussFilter.hpp
@@ -0,0 +1,81 @@
+//
+// Created by jstark on 2019-11-11.
+//
+/**
+ * @file GaussFilter.hpp
+ *
+ * @brief Header file containing function which applies Gaussian smoothing (Gaussian blur) to an n-dim. grid.
+ *
+ * @details In cases of very steep gradients of the initial (pre-redistancing) phi at the interface, it can be
+ * beneficial to smoothen phi prior to redistancing. If the initial function is a heaviside-like function, Gaussian
+ * smoothing transforms it into a rather sigmoid-like functions.
+ *
+ * @author Justina Stark
+ * @date November 2019
+ */
+#ifndef REDISTANCING_SUSSMAN_GAUSSFILTER_HPP
+#define REDISTANCING_SUSSMAN_GAUSSFILTER_HPP
+
+
+#include <iostream>
+#include <typeinfo>
+#include <cmath>
+
+#include "Grid/grid_dist_id.hpp"
+#include "HelpFunctions.hpp"
+#include "HelpFunctionsForGrid.hpp"
+
+/**@brief Discrete Gaussian filter with sigma=1.
+ *
+ * @details If Gaussian blur of sigma>1 is desired, just repeat smoothing.
+ *
+ * @tparam Phi_0_in Index of property that contains the values which should be Gaussian blurred.
+ * @tparam Phi_smoothed_out Index of property where the output of smoothing should be stored.
+ * @tparam grid_in_type Type of input grid.
+ * @param grid_in Input-OpenFPM-grid on which Gaussian blur should be applied.
+ */
+template <size_t Phi_0_in, size_t Phi_smoothed_out, typename grid_in_type>
+void gauss_filter(grid_in_type & grid_in)
+{
+ grid_in.template ghost_get<Phi_0_in>();
+
+ // create a temporary grid with same size on which smoothing is perfomed to not override Phi_0 of the input grid
+ typedef aggregate<double, double> props_temp;
+ typedef grid_dist_id<grid_in_type::dims, typename grid_in_type::stype, props_temp> g_temp_type;
+ g_temp_type g_temp(grid_in.getDecomposition(), grid_in.getGridInfoVoid().getSize(), Ghost<grid_in_type::dims, long int>(3));
+ // Some indices for better readability
+ static constexpr size_t Phi_0_temp = 0;
+ static constexpr size_t Phi_smoothed_temp = 1;
+
+ copy_gridTogrid<Phi_0_in, Phi_0_temp>(grid_in, g_temp, true); // copy incl. ghost
+
+ double gauss1d[] = {0.006, 0.061, 0.242, 0.383}; // 1D Gauss kernel for sigma = 1.0
+
+// loop over grid points. Gaussian kernel separated in its dimensions and first applied in x, then in y and z
+ for (int d = 0; d < g_temp_type::dims; d++)
+ {
+ g_temp.template ghost_get<Phi_0_temp, Phi_smoothed_temp>();
+ auto dom1 = g_temp.getDomainIterator();
+ while (dom1.isNext())
+ {
+ auto key = dom1.get();
+ auto key_g = g_temp.getGKey(key);
+
+ g_temp.template get<Phi_smoothed_temp>(key) =
+ gauss1d[0] * g_temp.template get<Phi_0_temp>(key.move(d, -3)) +
+ gauss1d[1] * g_temp.template get<Phi_0_temp>(key.move(d, -2)) +
+ gauss1d[2] * g_temp.template get<Phi_0_temp>(key.move(d, -1)) +
+ gauss1d[3] * g_temp.template get<Phi_0_temp>(key) +
+ gauss1d[2] * g_temp.template get<Phi_0_temp>(key.move(d, +1)) +
+ gauss1d[1] * g_temp.template get<Phi_0_temp>(key.move(d, +2)) +
+ gauss1d[0] * g_temp.template get<Phi_0_temp>(key.move(d, +3));
+ ++dom1;
+ }
+ copy_gridTogrid<Phi_smoothed_temp, Phi_0_temp>(g_temp, g_temp);
+ g_temp.template ghost_get<Phi_0_temp, Phi_smoothed_temp>();
+ }
+ copy_gridTogrid<Phi_smoothed_temp, Phi_smoothed_out>(g_temp, grid_in, true); // copy incl. ghost
+}
+
+
+#endif //REDISTANCING_SUSSMAN_GAUSSFILTER_HPP
diff --git a/src/level_set/redistancing_Sussman/HelpFunctions.hpp b/src/level_set/redistancing_Sussman/HelpFunctions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cd7d36730ecd32d94fc271efee3991340bda5f91
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/HelpFunctions.hpp
@@ -0,0 +1,84 @@
+//
+// Created by jstark on 2019-11-08.
+//
+/**
+ * @file HelpFunctions.hpp
+ *
+ * @brief Header file containing small mathematical help-functions that are needed e.g. for the Sussman redistancing.
+ *
+ * @author Justina Stark
+ * @date November 2019
+ */
+#ifndef REDISTANCING_SUSSMAN_HELPFUNCTIONS_HPP
+#define REDISTANCING_SUSSMAN_HELPFUNCTIONS_HPP
+
+#include <iostream>
+/**@brief Gets the sign of a variable.
+ *
+ * @tparam T Inferred type of input variable.
+ * @param val Scalar variable of arbitrary type for which the sign should be determined.
+ * @return Integer variable that contains a -1, if argument negative, 0 if argument=0, and +1 if argument positive.
+ */
+template <class T>
+int sgn(T val)
+{
+ return (T(0) < val) - (val < T(0));
+}
+/**@brief Gets the smoothed sign of a variable.
+ *
+ * @details Sign function with smoothing effect on numerical solution (see: <a href="https://www.math.uci
+ * .edu/~zhao/publication/mypapers/ps/local_levelset.ps">Peng \a et \a al. "A PDE-Based Fast Local
+ * Level Set Method"</a>, equation (36). Peng \a et \a al.: "The choice of approximation to S(d) by (36) solves the
+ * problem of the changing of sign of Phi (thus moving the interface across the cell boundary) in the reinitialization
+ * step when Phi is steep and speeds up the convergence when Phi is flat at the interface."
+ *
+ * @f[ sgn_{\epsilon}(\phi) = \frac{\phi}{ \sqrt{\phi^2 + |\nabla\phi|^2 \Delta x^2} }@f]
+ *
+ * where the \a φ is the approximation of the SDF from the last iteration.
+ *
+ * @tparam T Inferred type of the variable for which the smooth sign should be determined.
+ * @param val Scalar variable for which the smooth sign should be determined.
+ * @param epsilon Scalar variable containing the smoothing factor. In case of Peng: @f[\epsilon =
+ * |\nabla\phi|^2 \Delta x^2@f]
+ * @return Smooth sign of the argument variable: @f[sgn_{\epsilon}(\phi)@f].
+ */
+template <typename T>
+T smooth_S(T val, T epsilon)
+{
+ return (val / sqrt(val * val + epsilon * epsilon));
+}
+
+
+/**@brief Checks, if two values are sufficiently close to each other within a given tolerance range, as to be
+ * considered as approximately equal.
+ *
+ * @tparam T Inferred type of the two variables, for which it should be checked, if they are sufficiently close.
+ * @param val1 Variable that contains the first value.
+ * @param val2 Variable that contains the second value.
+ * @param tolerance Tolerance (or error) by which values are allowed to deviate while still being considered as equal.
+ * @return True, if \p val1 and \p val2 are the same +/- \p tolerance. False, if they differ by more
+ * than \p tolerance.
+ */
+template <class T>
+bool isApproxEqual(T val1, T val2, T tolerance)
+{
+ return (val1 <= val2 + tolerance && val1 >= val2 - tolerance);
+}
+
+/**@brief Appends the value of a given variable of any type to a textfile as string.
+ *
+ * @tparam T Inferred type of the input variable, whose value should be appended to textfile as string.
+ * @param textfile Std::string that contains the path and filename of the textfile, to which \p value should be added.
+ * @param value Variable that contains the value which should be appended to the \p textfile.
+ */
+template <typename T>
+void append_value_to_textfile(std::string & textfile, T value)
+{
+ std::ofstream out(textfile);
+ out << value;
+}
+
+
+
+
+#endif //REDISTANCING_SUSSMAN_HELPFUNCTIONS_HPP
diff --git a/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ed4a02bee2e7040e27b7ec8ae51990a12f2e6cb5
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
@@ -0,0 +1,231 @@
+//
+// Created by jstark on 2020-05-12.
+//
+/**
+ * @file HelpFunctionsForGrid.hpp
+ *
+ * @brief Header file containing help-functions that perform on OpenFPM-grids.
+ *
+ * @author Justina Stark
+ * @date May 2020
+ */
+#ifndef REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP
+#define REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP
+
+#include <iostream>
+#include <limits>
+
+#include "HelpFunctions.hpp"
+
+/**@brief Computes the time step for the iterative re-distancing fulfilling CFL condition.
+ *
+ * @tparam grid_type Inferred type of the input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Time step.
+ */
+template <typename grid_type>
+double get_time_step_CFL(grid_type &grid)
+{
+ double divisor = 0;
+ for (size_t d = 0; d < grid_type::dims; d++)
+ {
+ divisor += 1.0 / (grid.spacing(d) * grid.spacing(d));
+ }
+ return 1.0 / (2.0 * divisor);
+}
+
+/**@brief Initializes given property \p Prop of an OpenFPM grid including ghost layer with a given value from \p
+ * init_value.
+ *
+ * @tparam Prop Index of property that should be initialized with the value from \p init_value.
+ * @tparam grid_type Inferred type of input OpenFPM grid.
+ * @tparam T Inferred type of the variable containing the initialization value \p init_value.
+ * @param grid OpenFPM grid whose property \p Prop should be initialized, including its ghost layer.
+ * @param init_value Variable that contains the value that should be copied to \p Prop.
+ */
+template <size_t Prop, typename grid_type, typename T>
+void init_grid_and_ghost(grid_type & grid, T init_value)
+{
+ auto dom_ghost = grid.getDomainGhostIterator();
+ while(dom_ghost.isNext())
+ {
+ auto key = dom_ghost.get();
+ grid.template get<Prop>(key) = init_value;
+ ++dom_ghost;
+ }
+}
+
+/**@brief Copies the value stored in a given property from a given source grid to a given destination grid.
+ *
+ * @tparam attr_sc Index of property that contains the value that should be copied.
+ * @tparam attr_ds Index of property to which the value should be copied to.
+ * @tparam grid_source_type Inferred type of the source grid \p grid_sc.
+ * @tparam grid_dest_type Inferred type of the destination grid \p grid_ds.
+ * @param grid_sc OpenFPM grid from which value is copied (source).
+ * @param grid_ds OpenFPM grid to which value is copied to (destination).
+ * @param include_ghost Bool variable that defines if the copying should include the ghost layer. False, if not; true,
+ * if yes. Default: false.
+ */
+template <size_t attr_sc, size_t attr_ds, typename grid_source_type, typename grid_dest_type>
+void copy_gridTogrid(const grid_source_type & grid_sc, grid_dest_type & grid_ds, bool include_ghost=false)
+{
+ assert(grid_source_type::dims == grid_dest_type::dims);
+ assert(grid_sc.size() == grid_ds.size());
+ if (include_ghost)
+ {
+ auto dom_sc = grid_sc.getDomainGhostIterator();
+ auto dom_ds = grid_ds.getDomainGhostIterator();
+ while (dom_sc.isNext())
+ {
+ auto key_sc = dom_sc.get();
+ auto key_ds = dom_ds.get();
+ grid_ds.template get<attr_ds>(key_ds) = grid_sc.template get<attr_sc>(key_sc);
+ ++dom_sc;
+ ++dom_ds;
+ }
+ }
+ else
+ {
+ auto dom_sc = grid_sc.getDomainIterator();
+ auto dom_ds = grid_ds.getDomainIterator();
+ while (dom_sc.isNext())
+ {
+ auto key_sc = dom_sc.get();
+ auto key_ds = dom_ds.get();
+ grid_ds.template get<attr_ds>(key_ds) = grid_sc.template get<attr_sc>(key_sc);
+ ++dom_sc;
+ ++dom_ds;
+ }
+ }
+}
+
+/**@brief Determines the biggest spacing of a grid which is potentially anisotropic when comparing x, y (and z)
+ * coordinate axis.
+ *
+ * @tparam grid_type Inferred type of input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Double variable that contains the value of the biggest spacing.
+ */
+template <typename grid_type>
+double get_biggest_spacing(grid_type & grid)
+{
+ double h_max = 0;
+ for (size_t d = 0; d < grid_type::dims; d++)
+ {
+ if (grid.spacing(d) > h_max) h_max = grid.spacing(d);
+ }
+ return h_max;
+}
+
+/**@brief Determines the smallest spacing of a grid which is potentially anisotropic when comparing x, y (and z)
+ * coordinate axis.
+ *
+ * @tparam grid_type Inferred type of input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Double variable that contains the value of the smallest spacing.
+ */
+template <typename grid_type>
+double get_smallest_spacing(grid_type & grid)
+{
+ double spacing [grid_type::dims];
+ for (size_t d = 0; d < grid_type::dims; d++)
+ {
+ spacing[d] = grid.spacing(d);
+ }
+ std::sort(std::begin(spacing), std::end(spacing));
+ return spacing[0];
+}
+
+/**@brief Computes the average difference between the values stored at two different properties of the same grid, that
+ * is, the total difference of these values summed over all the grid nodes divided by the gridsize.
+ *
+ * @tparam Prop1 Index of the first property.
+ * @tparam Prop2 Index of the second property.
+ * @tparam grid_type Inferred type of the input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Double variable that contains the sum over all grid nodes of the difference between the value stored at \p
+ * Prop1 and the value stored at \p Prop2.
+ */
+template <size_t Prop1, size_t Prop2, typename grid_type>
+double average_difference(grid_type & grid)
+{
+ double total_diff = 0;
+ auto dom = grid.getDomainIterator();
+ while (dom.isNext())
+ {
+ auto key = dom.get();
+ total_diff += abs( grid.template get<Prop1>(key) - grid.template get<Prop2>(key) );
+ ++dom;
+ }
+ return total_diff / grid.size();
+}
+
+/**@brief Determines the maximum value stored on a given grid at a given property.
+ *
+ * @tparam Prop Index of property for which maximum should be found.
+ * @tparam grid_type Inferred type of the input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Double variable that contains the maximum value of property \p Prop in \p grid.
+ */
+template <size_t Prop, typename grid_type>
+double get_max_val(grid_type & grid)
+{
+ double max_value = std::numeric_limits<double>::lowest();
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ if (grid.template get<Prop>(key) > max_value) max_value = grid.template get<Prop>(key);
+ ++dom;
+ }
+ auto & v_cl = create_vcluster();
+ v_cl.max(max_value);
+ v_cl.execute();
+ return max_value;
+}
+
+/**@brief Determines the minimum value stored on a given grid at a given property.
+ *
+ * @tparam Prop Index of property for which minimum should be found.
+ * @tparam grid_type Inferred type of the input grid.
+ * @param grid Input OpenFPM grid.
+ * @return Double variable that contains the minimum value of property \p Prop in \p grid.
+ */
+template <size_t Prop, typename grid_type>
+double get_min_val(grid_type & grid)
+{
+ double min_value = std::numeric_limits<double>::max();
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ if (grid.template get<Prop>(key) < min_value) min_value = grid.template get<Prop>(key);
+ ++dom;
+ }
+ auto & v_cl = create_vcluster();
+ v_cl.min(min_value);
+ v_cl.execute();
+ return min_value;
+}
+
+/**@brief Determines the sign of a value stored at a given property and stores it in another property.
+ *
+ * @tparam Prop_in Index of property that contains the value whose sign should be determined.
+ * @tparam Prop_out Index of property where the sign should be written to.
+ * @tparam grid_type Inferred type of the input grid.
+ * @param grid Input OpenFPM grid.
+ */
+template <size_t Prop_in, size_t Prop_out, typename grid_type>
+void init_sign_prop(grid_type & grid)
+{
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ grid.template get<Prop_out>(key) = sgn(grid.template get<Prop_in>(key));
+ ++dom;
+ }
+}
+
+
+#endif //REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP
diff --git a/src/level_set/redistancing_Sussman/NarrowBand.hpp b/src/level_set/redistancing_Sussman/NarrowBand.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6532defc6b7b91a88ae271bd7af6f798fe34969
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/NarrowBand.hpp
@@ -0,0 +1,397 @@
+//
+// Created by jstark on 2020-05-16.
+//
+/**
+ * @file NarrowBand.hpp
+ * @class NarrowBand
+ *
+ * @brief Class for getting the narrow band around the interface.
+ *
+ * @details Places particle at those grid point positions, which lay inside a narrow band of user-defined width
+ * around the interface. The respective SDF or arbitrary other property value is then also copied from the grid node
+ * to that particle occupying the same position.
+ *
+ *
+ * @author Justina Stark
+ * @date May 2020
+ */
+
+#ifndef REDISTANCING_SUSSMAN_NARROWBAND_HPP
+#define REDISTANCING_SUSSMAN_NARROWBAND_HPP
+
+// Include standard library header files
+#include <iostream>
+
+// Include OpenFPM header files
+#include "Vector/vector_dist.hpp"
+#include "Grid/grid_dist_id.hpp"
+#include "data_type/aggregate.hpp"
+#include "Decomposition/CartDecomposition.hpp"
+
+// Include level-set-method related header files
+#include "ComputeGradient.hpp"
+
+/**@brief Class for getting the narrow band around the interface
+ * @file NarrowBand.hpp
+ * @class NarrowBand
+ * @tparam grid_in_type Inferred type of input grid that stores the signed distance function Phi_SDF.
+ */
+template <typename grid_in_type>
+class NarrowBand
+{
+public:
+ /** @brief Constructor taking the thickness of the narrow band as #grid points in order to initialize the
+ * lower and upper bound for the narrow band. Initializes the temporary internal grid.
+ *
+ * @param grid_in Input grid with min. 1 property: Phi_SDF.
+ * @param thickness Width of narrow band in # grid points.
+ */
+ NarrowBand(const grid_in_type & grid_in,
+ size_t thickness) // thickness in # grid points
+ : g_temp(grid_in.getDecomposition(), grid_in.getGridInfoVoid().getSize(), Ghost<grid_in_type::dims, long int>(3))
+ {
+ set_bounds(thickness, grid_in);
+ }
+ /** @brief Constructor taking the thickness of the narrow band as physical width in space in order to initialize the
+ * lower and upper bound for the narrow band. Initializes the temporary internal grid.
+ *
+ * @param grid_in Input grid with min. 1 property: Phi_SDF.
+ * @param thickness Width of narrow band as physical width.
+ */
+ NarrowBand(const grid_in_type & grid_in,
+ double thickness) // thickness as physical width
+ : g_temp(grid_in.getDecomposition(), grid_in.getGridInfoVoid().getSize(), Ghost<grid_in_type::dims, long int>(3))
+ {
+ set_bounds(thickness, grid_in);
+ }
+ /** @brief Constructor taking the inside and outside physical width of the narrow band in order to initialize the
+ * lower and upper bound for the narrow band. Initializes the temporary internal grid.
+ *
+ * @param grid_in Input grid with min. 1 property: Phi_SDF.
+ * @param width_outside Extension of narrow band as physical width inside of object (Phi > 0).
+ * @param width_inside Extension of narrow band as physical width outside of object (Phi < 0).
+ */
+ NarrowBand(const grid_in_type & grid_in,
+ double width_outside, // physical width of nb inside of object -> Phi > 0
+ double width_inside) // physical width nb outside of object -> Phi < 0
+ : g_temp(grid_in.getDecomposition(), grid_in.getGridInfoVoid().getSize(), Ghost<grid_in_type::dims, long int>(3))
+ {
+ set_bounds(width_outside, width_inside, grid_in);
+ }
+
+
+ /**@brief Aggregated properties for the temporary grid.
+ *
+ * @details The properties are Phi_SDF (result from redistancing), gradient of Phi, L2 norm of gradient of Phi
+ * (=gradient magnitude) and sign of Phi_SDF (for the upwinding).
+ * */
+ typedef aggregate<double, double[grid_in_type::dims], double, int> props_temp;
+ /** @brief Type definition for the temporary grid.
+ */
+ typedef grid_dist_id<grid_in_type::dims, typename grid_in_type::stype, props_temp> g_temp_type;
+ /**
+ * @brief Create temporary grid, which is only used inside the class to get the gradients.
+ *
+ * @details The grid is needed, when the user wants a narrow band which also contains the upwind gradient of phi.
+ * It contains the following 5 properties: Phi_SDF (result from redistancing), gradient of Phi, L2 norm of gradient
+ * of Phi (=gradient magnitude) and sign of Phi_SDF (for the upwinding).
+ * */
+ g_temp_type g_temp;
+
+ /** @brief Calls NarrowBand::get_narrow_band_phi_only which fills the empty vector passed as argument with particles
+ * by placing these within a narrow band around the interface. Only the SDF is copied from the grid properties to
+ * the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param vd Empty vector with same spatial scaling (box) as the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, typename vector_type, typename grid_type>
+ void get_narrow_band(grid_type & grid, vector_type & vd)
+ {
+ get_narrow_band_phi_only<Phi_SDF_grid, Phi_SDF_vd>(grid, vd);
+ }
+ /** @brief Calls NarrowBand::get_narrow_band_with_gradients which fills the empty vector passed as argument with particles by placing these within a narrow band around the
+ * interface. SDF and Phi_grad_temp are copied from the temp. grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam Phi_grad Index of property that should store the gradient of phi in the narrow band particle vector.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param vd Empty vector with same spatial scaling (box) as the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, size_t Phi_grad, typename vector_type, typename grid_type>
+ void get_narrow_band(grid_type & grid, vector_type & vd)
+ {
+ get_narrow_band_with_gradients<Phi_SDF_grid, Phi_SDF_vd, Phi_grad>(grid, vd);
+ }
+ /** @brief Calls NarrowBand::get_narrow_band_with_gradients which fills the empty vector passed as argument with
+ * particles by placing these within a narrow band around the interface. SDF, Phi_grad_temp and
+ * Phi_magnOfGrad_temp are copied from the temp. grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam Phi_grad Index of property that should store the gradient of phi in the narrow band particle vector.
+ * @tparam Phi_magnOfGrad Index of property that should store the gradient magnitude of phi in the narrow band particles.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param vd Empty vector with same spatial scaling (box) as the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, size_t Phi_grad, size_t Phi_magnOfGrad, typename vector_type, typename grid_type>
+ void get_narrow_band(grid_type & grid, vector_type & vd)
+ {
+ get_narrow_band_with_gradients<Phi_SDF_grid, Phi_SDF_vd, Phi_grad, Phi_magnOfGrad>(grid, vd);
+ }
+ /** @brief Calls NarrowBand::get_narrow_band_one_prop which fills the empty vector passed as argument with
+ * particles by placing these within a narrow band around the interface. An arbitrary property is copied from the
+ * grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Prop1_grid Index of arbitrary grid property that should be copied to the particles (prop. value must
+ * be a scalar).
+ * @tparam Prop1_vd Index of particle property, where grid property should be copied to (prop. value must be a
+ * scalar).
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param grid Grid of arb. dims. storing the SDF (result of redistancing) and any arbitrary other (scalar)
+ * property.
+ * @param vd Empty vector with same spatial scaling (box) as the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Prop1_grid, size_t Prop1_vd, typename vector_type, typename grid_type>
+ void get_narrow_band_copy_specific_property(grid_type & grid, vector_type & vd)
+ {
+ get_narrow_band_one_prop<Phi_SDF_grid, Prop1_grid, Prop1_vd>(grid, vd);
+ }
+
+private:
+ // Some indices for better readability
+ static const size_t Phi_SDF_temp = 0; ///< Property index of Phi_SDF on the temporary grid.
+ static const size_t Phi_grad_temp = 1; ///< Property index of gradient of Phi on the temporary grid.
+ static const size_t Phi_magnOfGrad_temp = 2; ///< Property index of gradient magnitude of Phi on the temp. grid.
+ static const size_t Phi_sign_temp = 3; ///< Property index of sign of Phi on the temporary grid.
+
+ double b_low; ///< Narrow band extension towards the object outside.
+ double b_up; ///< Narrow band extension towards the object inside.
+
+ /** @brief Set the member variable NarrowBand::b_low and NarrowBand::b_up.
+ *
+ * @param thickness Thickness of narrow band in number of grid points.
+ * @param grid_in Input grid storing the signed distance function Phi_SDF (redistancing output).
+ */
+ void set_bounds(size_t thickness, const grid_in_type & grid_in)
+ {
+ b_low = - ceil(thickness / 2.0) * get_biggest_spacing(grid_in);
+ b_up = ceil(thickness / 2.0) * get_biggest_spacing(grid_in);
+ }
+ /** @brief Set the member variable NarrowBand::b_low and NarrowBand::b_up.
+ *
+ * @param thickness Thickness of narrow band as physical width.
+ * @param grid_in Input grid storing the signed distance function Phi_SDF (redistancing output).
+ */
+ void set_bounds(double thickness, const grid_in_type & grid_in)
+ {
+ b_low = - thickness / 2.0;
+ b_up = thickness / 2.0;
+ }
+ /** @brief Set the member variable NarrowBand::b_low and NarrowBand::b_up.
+ *
+ * @param lower_bound Extension of narrow band as physical width inside of object (Phi > 0).
+ * @param upper_bound Extension of narrow band as physical width outside of object (Phi < 0).
+ * @param grid_in Input grid storing the signed distance function Phi_SDF (redistancing output).
+ */
+ void set_bounds(double lower_bound, double upper_bound, const grid_in_type & grid_in)
+ {
+ b_low = lower_bound;
+ b_up = upper_bound;
+ }
+
+ /**@brief Initialize the internal temporary grid.
+ *
+ * @details Copies Phi_SDF from the input grid to the temorary grid.
+ *
+ * @tparam Phi_SDF Index of property storing the signed distance function in the input grid.
+ * @param grid_in Input grid storing the signed distance function Phi_SDF (redistancing output).
+ */
+ template <size_t Phi_SDF>
+ void initialize_temporary_grid(const grid_in_type & grid_in)
+ {
+ copy_gridTogrid<Phi_SDF, Phi_SDF_temp>(grid_in, g_temp); // Copy Phi_SDF from the input grid to the temorary grid
+ init_sign_prop<Phi_SDF_temp, Phi_sign_temp>(g_temp); // initialize Phi_sign_temp with the sign of the
+ // input Phi_SDF
+ get_upwind_gradient<Phi_SDF_temp, Phi_sign_temp, Phi_grad_temp>(g_temp); // Get initial gradients
+ get_gradient_magnitude<Phi_grad_temp, Phi_magnOfGrad_temp>(g_temp); // Get initial magnitude of gradients
+ }
+ /**@brief Checks if a value for Phi_SDF lays within the narrow band.
+ *
+ * @param Phi Value of the signed distance function Phi_SDF.
+ * @return True, if point lays within the narrow band; False if not.
+ */
+ bool within_narrow_band(double Phi)
+ {
+ return (Phi >= b_low && Phi <= b_up);
+ }
+ /** @brief Fills the empty vector passed as argument with particles by placing these within a narrow band around the
+ * interface. Only the SDF is copied from the grid properties to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param[in] grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param[in] vd Empty vector with same spatial scaling (box) as the grid.
+ * @param[out] vd Particle vector containing the narrow-band-particles that store the SDF provided by the grid.
+ *
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, typename grid_type, typename vector_type>
+ void get_narrow_band_phi_only(grid_type & grid, vector_type & vd)
+ {
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_g = grid.getGKey(key);
+ if(within_narrow_band(grid.template get<Phi_SDF_grid>(key)))
+ {
+ // add particle to vd_narrow_band
+ vd.add();
+ // assign coordinates and properties from respective grid point to particle
+ for(size_t d = 0; d < grid_type::dims; d++)
+ {
+ vd.getLastPos()[d] = key_g.get(d) * grid.getSpacing()[d];
+ vd.template getLastProp<Phi_SDF_vd>() = grid.template get<Phi_SDF_grid>(key);
+ }
+ }
+ ++dom;
+ }
+ }
+ /** @brief Fills the empty vector passed as argument with particles by placing these within a narrow band around the
+ * interface. SDF and Phi_grad_temp are copied from the temp. grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam Phi_grad Index of property that should store the gradient of phi in the narrow band particle vector.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param[in] grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param[in] vd Empty vector with same spatial scaling (box) as the grid.
+ * @param[out] vd Particle vector containing the narrow-band-particles that store the SDF and Phi_grad provided
+ * by the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, size_t Phi_grad, typename grid_type, typename vector_type>
+ void get_narrow_band_with_gradients(grid_type & grid, vector_type & vd)
+ {
+ initialize_temporary_grid<Phi_SDF_grid>(grid);
+ auto dom = g_temp.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_g = g_temp.getGKey(key);
+ if(within_narrow_band(g_temp.template get<Phi_SDF_temp>(key)))
+ {
+ // add particle to vd_narrow_band
+ vd.add();
+ // assign coordinates and properties from respective grid point to particle
+ for(size_t d = 0; d < grid_type::dims; d++)
+ {
+ vd.getLastPos()[d] = key_g.get(d) * g_temp.getSpacing()[d];
+ vd.template getLastProp<Phi_SDF_vd>() = g_temp.template get<Phi_SDF_temp>(key);
+ vd.template getLastProp<Phi_grad>()[d] = g_temp.template get<Phi_grad_temp>(key)[d];
+ }
+ }
+ ++dom;
+ }
+ }
+ /** @brief Fills the empty vector passed as argument with particles by placing these within a narrow band around the
+ * interface. SDF, Phi_grad_temp and Phi_magnOfGrad_temp are copied from the temp. grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Phi_SDF_vd Index of property that should store the SDF in the narrow band particle vector.
+ * @tparam Phi_grad Index of property that should store the gradient of phi in the narrow band particle vector.
+ * @tparam Phi_magnOfGrad Index of property that should store the gradient magnitude of phi in the narrow band particles.
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param[in] grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param[in] vd Empty vector with same spatial scaling (box) as the grid.
+ * @param[out] vd Particle vector containing the narrow-band-particles that store the SDF, Phi_grad and
+ * Phi_magnOfGrad provided by the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Phi_SDF_vd, size_t Phi_grad, size_t Phi_magnOfGrad, typename grid_type, typename vector_type>
+ void get_narrow_band_with_gradients(grid_type & grid, vector_type & vd)
+ {
+ initialize_temporary_grid<Phi_SDF_grid>(grid);
+ auto dom = g_temp.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_g = g_temp.getGKey(key);
+ if(within_narrow_band(g_temp.template get<Phi_SDF_temp>(key)))
+ {
+ // add particle to vd_narrow_band
+ vd.add();
+ // assign coordinates and properties from respective grid point to particle
+ for(size_t d = 0; d < grid_type::dims; d++)
+ {
+ vd.getLastPos()[d] = key_g.get(d) * g_temp.getSpacing()[d];
+ vd.template getLastProp<Phi_SDF_vd>() = g_temp.template get<Phi_SDF_temp>(key);
+ vd.template getLastProp<Phi_grad>()[d] = g_temp.template get<Phi_grad_temp>(key)[d];
+ vd.template getLastProp<Phi_magnOfGrad>() = g_temp.template get<Phi_magnOfGrad_temp>(key);
+ }
+ }
+ ++dom;
+ }
+ }
+ /** @brief Fills the empty vector passed as argument with particles by placing these within a narrow band around the
+ * interface. An arbitrary property is copied from the grid to the respective particle.
+ *
+ * @tparam Phi_SDF_grid Index of property storing the signed distance function in the input grid.
+ * @tparam Prop1_grid Index of arbitrary grid property that should be copied to the particles (prop. value must
+ * be a scalar).
+ * @tparam Prop1_vd Index of particle property, where grid property should be copied to (prop. value must be a
+ * scalar).
+ * @tparam vector_type Inferred type of the particle vector.
+ * @tparam grid_type Inferred type of the grid storing the SDF.
+ *
+ * @param[in] grid Grid of arb. dims. storing the SDF (result of redistancing).
+ * @param[in] vd Empty vector with same spatial scaling (box) as the grid.
+ * @param[out] vd Particle vector containing the narrow-band-particles that store an arbitrary property provided by
+ * the grid.
+ */
+ template <size_t Phi_SDF_grid, size_t Prop1_grid, size_t Prop1_vd, typename grid_type, typename vector_type>
+ void get_narrow_band_one_prop(grid_type & grid, vector_type & vd)
+ {
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_g = grid.getGKey(key);
+ if(within_narrow_band(grid.template get<Phi_SDF_grid>(key)))
+ {
+ // add particle to vd_narrow_band
+ vd.add();
+ // assign coordinates and properties from respective grid point to particle
+ for(size_t d = 0; d < grid_type::dims; d++)
+ {
+ vd.getLastPos()[d] = key_g.get(d) * grid.getSpacing()[d];
+ vd.template getLastProp<Prop1_vd>() = grid.template get<Prop1_grid>(key);
+ }
+ }
+ ++dom;
+ }
+ }
+};
+
+
+#endif //REDISTANCING_SUSSMAN_NARROWBAND_HPP
diff --git a/src/level_set/redistancing_Sussman/RedistancingSussman.hpp b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..397e84013dcd6474477c1fbcf02522bee30159a5
--- /dev/null
+++ b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
@@ -0,0 +1,545 @@
+//
+// Created by jstark on 2020-04-06.
+//
+/**
+ * @file RedistancingSussman.hpp
+ * @class RedistancingSussman
+ *
+ * @brief Class for reinitializing a level-set function into a signed distance function using Sussman redistancing.
+ *
+ * @details First-time-only redistancing step (see: <a href="https://www.researchgate
+ * .net/publication/2642502_An_Efficient_Interface_Preserving_Level_Set_Re
+ * -Distancing_Algorithm_And_Its_Application_To_Interfacial_Incompressible_Fluid_Flow.html">M. Sussman and E. Fatemi,
+ * “Efficient, interface-preserving level set redistancing algorithm and its application to interfacial
+ * incompressible fluid flow†(1999)</a> ). In the Sussman-redistancing, a level-set-function Phi_0, which can be a
+ * simple step-function having positive values inside the object and negative values outside, is reinitialized to a
+ * signed distance function Phi_SDF by finding the steady-state solution of the following PDE:
+ * @f[ \phi_{t} + sgn_{\epsilon}(\phi)(|\nabla\phi| - 1) = 0 @f]
+ *
+ * The signed distance function is defined as:
+ *
+ * @f[ \phi_{\text{SDF}} = \begin{cases}
+ * +d & \text{orthogonal distance to the closest surface of a point lying inside the object} \\
+ * 0 & \text{object surface} \\
+ * -d & \text{orthogonal distance to the closest surface of a point lying outside the object} \\
+ * \end{cases} @f]
+ *
+ * @author Justina Stark
+ * @date April 2020
+ */
+
+#ifndef REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP
+#define REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP
+
+// Include standard library header files
+#include <iostream>
+#include <string>
+#include <fstream>
+// Include OpenFPM header files
+#include "Vector/vector_dist.hpp"
+#include "Grid/grid_dist_id.hpp"
+#include "data_type/aggregate.hpp"
+#include "Decomposition/CartDecomposition.hpp"
+
+// Include self-written header files
+#include "HelpFunctions.hpp"
+#include "HelpFunctionsForGrid.hpp"
+#include "GaussFilter.hpp"
+#include "ComputeGradient.hpp"
+
+/** @brief Optional convergence criterium checking the total change.
+ *
+ * @details The change between the current and the previous iterations is computed as sum over all the If
+ * Conv_tol_change
+ * .check = true, the total change of Phi
+ * between the iterations is considered and the
+ * redistancing finishes when the Conv_tol_change.value (or the max-iteration) is reached.
+ */
+struct Conv_tol_change
+{
+ bool check = true; ///< If true, the total change of Phi (see DistFromSol::change)
+ ///< between the iterations is considered and the redistancing finishes when the Conv_tol_change.value (or the
+ ///< max-iteration) is reached. If false, the change is not considered as steady-state
+ ///< criterion. Default: true.
+ double value = 1e-6; ///< Double variable that stores the tolerance value for the change at which the iterative
+ ///< redistancing process is considered as steady-state. (see also DistFromSol::change)
+};
+
+/** @brief Optional convergence criterium checking the residual.
+ *
+ * @details If Conv_tol_residual.check = true, the residual, that is abs(magnitude gradient of phi - 1), is considered
+ * and the redistancing finishes when the Conv_tol_residual.value (or the max-iteration) is reached.
+ */
+struct Conv_tol_residual
+{
+ bool check = false; ///< If true, the residual of Phi (see DistFromSol::residual) is considered and the
+ ///< redistancing finishes when the Conv_tol_residual.value (or the
+ ///< max-iteration) is reached. If false, the change is not considered as steady-state criterion. Default: false.
+ double value = 1e-1; ///< Double variable that stores the tolerance value for the residual at which the
+ ///< iterative redistancing process is considered as steady-state. (see also
+ ///< DistFromSol::residual)
+};
+
+/** @brief Structure to bundle options for redistancing.
+ * @struct Redist_options
+ * @details For the redistancing, we can choose some options. These options will then be passed bundled as a structure to
+ * the redistancing function. Setting these options is optional, since they all have a Default value as well.
+ *
+ * @param sigma: Sigma of the gaussian kernel, which is used for gaussian smooting Phi_0. If the
+ * initial gradient of
+ * phi_0 at the interface is too large and no sigma is chosen or chosen too small, gauss smoothing will
+ * automatically be applied until phi gradient magnitude <= 12, regardless of which sigma is chosen by
+ * the user. Default = 0.
+ * @param min_iter: Minimum number of iterations before steady state in narrow band will be checked (Default: 100).
+ * @param max_iter: Maximum number of iterations you want to run the redistancing, even if steady state might not yet
+ * have been reached (Default: 1e6).
+ * @param convTolChange.value: Convolution tolerance for the normalized total change of Phi in the narrow band between
+ * two consecutive iterations (Default: 1e-6).
+ * @param convTolChange.check: Set true, if you want to use the normalized total change between two iterations as
+ * measure of how close you are to the steady state solution. Redistancing will then stop
+ * if convTolChange.value is reached or if the current iteration is bigger than max_iter.
+ * @param convTolResidual.value: Convolution tolerance for the residual, that is abs(magnitude gradient of phi - 1) of
+ * Phi in the narrow band (Default 1e-1).
+ * @param convTolResidual.check: Set true, if you want to use the residual of the current iteration as measure of how
+ * close you are to the steady state solution. Redistancing will then stop if
+ * convTolResidual.value is reached or if the current iteration is bigger than max_iter.
+ * @param interval_check_convergence: Interval of number of iterations at which convergence to steady state is checked
+ * (Default: 100).
+ * @param width_NB_in_grid_points: Width of narrow band in number of grid points. Must be at least 4, in order to
+ * have at least 2 grid points on each side of the interface. Is automatically set
+ * to 4, if a value smaller than 4 is chosen (Default: 4).
+ * @param print_current_iterChangeResidual: If true, the number of the current iteration, the corresponding change
+ * w.r.t the previous iteration and the residual is printed (Default: false).
+ * @param print_steadyState_iter: If true, the number of the steady-state-iteration, the corresponding change
+ * w.r.t the previous iteration and the residual is printed (Default: false).
+ */
+struct Redist_options
+{
+ size_t sigma = 0; // if you want to apply Gaussian smoothing prior to the redistancing. Smoothing is applied
+ // automatically if h_min <= 1e-6 (otherwise gradient at interface too big)
+ size_t min_iter = 1000;
+ size_t max_iter = 1e6;
+
+ Conv_tol_change convTolChange;
+ Conv_tol_residual convTolResidual;
+
+ size_t interval_check_convergence = 100;
+ size_t width_NB_in_grid_points = 4;
+ bool print_current_iterChangeResidual = false;
+ bool print_steadyState_iter = true;
+};
+
+/** @brief Bundles total residual and total change over all the grid points.
+ *
+ * @details In order to evaluate, how close the current solution is away from the convergence criterium.
+ *
+ * @see RedistancingSussman::get_residual_and_change_NB()
+ */
+struct DistFromSol
+{
+ double residual; ///< Double variable that contains the absolute value of how far the gradient magnitude of
+ ///< the current \a φ of iteration number \a i is away from being equal to 1: @f[ abs(|\nabla\phi_i| - 1 ) @f]
+ ///< It is computed for all grid points that lie within the narrow band and
+ ///< normalized to the number of grid points that lie in that narrow band.
+ double change; ///< Double variable that contains the absolute value of the change of \a φ between the
+ ///< current
+ ///< iteration \a i and the previous iteration \a i-1: @f[ abs( \phi_i - \phi_{i-1} ) @f] It is
+ ///< computed for all grid
+ ///< points that lie within the narrow band and normalized to the number of grid points that
+ ///< lie in that narrow band.
+};
+
+
+/**@brief Class for reinitializing a level-set function into a signed distance function using Sussman redistancing.
+ * @file RedistancingSussman.hpp
+ * @class RedistancingSussman
+ * @tparam grid_in_type Inferred type of input grid, which stores the initial level-set function Phi_0.
+ */
+template <typename grid_in_type>
+class RedistancingSussman
+{
+public:
+ /** @brief Constructor initializing the redistancing options, the temporary internal grid and reference variable
+ * to the input grid.
+ *
+ * @param grid_in Input grid with min. 2 properties: 1.) Phi_0, 2.) Phi_SDF <- will be overwritten with
+ * re-distancing result
+ * @param redistOptions User defined options for the Sussman redistancing process
+ *
+ */
+ RedistancingSussman(grid_in_type &grid_in, Redist_options &redistOptions) : redistOptions(redistOptions),
+ r_grid_in(grid_in),
+ g_temp(grid_in.getDecomposition(),
+ grid_in.getGridInfoVoid().getSize(),
+ Ghost<grid_in_type::dims, long int>(
+ 3))
+ {
+ time_step = get_time_step_CFL(g_temp);
+ assure_minimal_thickness_of_NB();
+ }
+
+ /**@brief Aggregated properties for the temporary grid.
+ *
+ * @details The initial (input) Phi_0 (will be updated by Phi_{n+1} after each redistancing step),
+ * Phi_{n+1} (received from redistancing),
+ * gradient of Phi_{n+1},
+ * L2 norm of gradient of Phi_{n+1} (=gradient magnitude),
+ * sign of the original input Phi_0 (for the upwinding).
+ */
+ typedef aggregate<double, double, double[grid_in_type::dims], double, int> props_temp;
+ /** @brief Type definition for the temporary grid.
+ */
+ typedef grid_dist_id<grid_in_type::dims, typename grid_in_type::stype, props_temp> g_temp_type;
+ /**
+ * @brief Create temporary grid, which is only used inside the class for the redistancing.
+ *
+ * @details The temporary grid stores the following 5 properties:
+ * the initial (input) Phi_0 (will be updated by Phi_{n+1} after each redistancing step),
+ * Phi_{n+1}(received from redistancing),
+ * gradient of Phi_{n+1},
+ * L2 norm of gradient of Phi_{n+1} (=gradient magnitude),
+ * sign of the original input Phi_0 (for the upwinding).
+ */
+ g_temp_type g_temp;
+
+ /**@brief Runs the Sussman-redistancing.
+ *
+ * @details Copies Phi_0 from input grid to an internal temporary grid which allows
+ * having more properties. Computes the gradients. Applies gauss smoothing if sigma set in redistOptions or if
+ * initial gradients to steep. Runs the redistancing on the internal tenporary grid. Copies resulting signed
+ * distance function to the Phi_SDF_out property of the input grid.
+ */
+ template<size_t Phi_0_in, size_t Phi_SDF_out> void run_redistancing()
+ {
+ init_temp_grid<Phi_0_in>();
+ if (redistOptions.sigma > 0)
+ {
+ run_gauss_filter(redistOptions.sigma);
+ }
+ init_sign_prop<Phi_0_temp, Phi_0_sign_temp>(
+ g_temp); // initialize Phi_0_sign_temp with the sign of the initial (pre-redistancing) Phi_0
+ get_upwind_gradient<Phi_0_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp); // Get initial gradients
+ get_gradient_magnitude<Phi_grad_temp, Phi_magnOfGrad_temp>(g_temp); // Get initial magnitude of gradients
+
+ smoothing_if_grad_steep();
+
+ iterative_redistancing(
+ g_temp); // Do the redistancing on the temporary grid
+ copy_gridTogrid<Phi_nplus1_temp, Phi_SDF_out>(g_temp, r_grid_in); // Copy resulting SDF function to
+ // input grid
+ }
+
+ /** @brief Overwrite the time_step found via CFL condition with an individual time_step.
+ *
+ * @details If the user wants to overwrite the time_step found via CFL condition with an individual time_step.
+ * Should be only used carefully, time_step must not be too large (jump over solution) nor too small (extremely
+ * slow).
+ *
+ * @param dt Artificial time step by which re-distancing should be performed.
+ *
+ */
+ template<typename T>
+ void set_user_time_step(T dt)
+ {
+ time_step = dt;
+ }
+ /** @brief Access the artificial timestep (private member) which will be used for the iterative redistancing.
+ * @see get_time_step_CFL(g_temp_type &grid), set_user_time_step()
+ */
+ double get_time_step()
+ {
+ /// This timestep is computed according to the grid spacing fulfilling the CFL condition.
+ return time_step;
+ }
+
+private:
+ // Some indices for better readability
+ static constexpr size_t Phi_0_temp = 0; ///< Property index of Phi_0 on the temporary grid.
+ static constexpr size_t Phi_nplus1_temp = 1; ///< Property index of Phi_n+1 on the temporary grid.
+ static constexpr size_t Phi_grad_temp = 2; ///< Property index of gradient of Phi_n on the temporary grid.
+ static constexpr size_t Phi_magnOfGrad_temp = 3; ///< Property index of gradient magn. of Phi_n (temporary grid).
+ static constexpr size_t Phi_0_sign_temp = 4; ///< Property index of sign of initial (input) Phi_0 (temp. grid).
+
+
+ // Member variables
+ Redist_options redistOptions; ///< Instantiate redistancing options.
+ grid_in_type &r_grid_in; ///< Define reference to input grid.
+
+ double h_max = get_biggest_spacing(g_temp); ///< Grid spacing in less resolved direction.
+ double h_min = get_smallest_spacing(g_temp); ///< Grid spacing in higher resolved direction.
+ /// Transform the half-bandwidth in no_of_grid_points into physical half-bandwidth kappa.
+ double kappa = ceil(redistOptions.width_NB_in_grid_points / 2.0) * h_max;
+ /**@brief Artificial timestep for the redistancing iterations.
+ * @see get_time_step_CFL(g_temp_type &grid), get_time_step(), set_user_time_step()
+ */
+ double time_step;
+
+ // Member functions
+ /** @brief Copies values from input grid to internal temporary grid and initializes ghost layer with minimum
+ * value of input grid.
+ */
+ template<size_t Phi_0_in>
+ void init_temp_grid()
+ {
+ double min_value = get_min_val<Phi_0_in>(r_grid_in); // get minimum Phi_0 value on the input grid
+ init_grid_and_ghost<Phi_0_temp>(g_temp, min_value); // init. Phi_0_temp (incl. ghost) with min. Phi_0
+ init_grid_and_ghost<Phi_nplus1_temp>(g_temp, min_value); // init. Phi_nplus1_temp (incl. ghost) with min. Phi_0
+ copy_gridTogrid<Phi_0_in, Phi_0_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_0_temp
+ }
+
+ /** @brief Checks if narrow band thickness >= 4 grid points. Else, sets it to 4 grid points.
+ *
+ * @details Makes sure, that the narrow band within which the convergence criteria are checked during the
+ * redistancing, is thick enough.
+ */
+ void assure_minimal_thickness_of_NB()
+ {
+ if (redistOptions.width_NB_in_grid_points < 4)
+ {
+ redistOptions.width_NB_in_grid_points = 4;
+ } // overwrite kappa if set too small by user
+ }
+
+ /** @brief Convolves with a Gauss kernel if initial gradient at the interface too steep.
+ */
+ void smoothing_if_grad_steep()
+ {
+ double max_grad = get_max_val<Phi_magnOfGrad_temp>(g_temp);
+ while (max_grad > 12)
+ {
+ run_gauss_filter(1);
+ // Update gradient
+ get_upwind_gradient<Phi_0_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp);
+ get_gradient_magnitude<Phi_grad_temp, Phi_magnOfGrad_temp>(g_temp);
+ // Update steepest gradient
+ max_grad = get_max_val<Phi_magnOfGrad_temp>(g_temp);
+ }
+ }
+ /** @brief Finds number of iterations a gauss kernel of sigma=1 has to be applied for getting an equivalent to
+ * convolution with a sigma=input sigma.
+ *
+ * @param sigma Sigma of gauss kernel with which grid should be smoothed prior to re-distancing.
+ *
+ * @return Number of iterations required with gauss of sigma = 1.
+ *
+ */
+ inline size_t count_gauss_repitition(const size_t sigma)
+ {
+ return sigma * sigma;
+ }
+
+ /** @brief Convolves the grid with a gauss filter.
+ *
+ * @param sigma Sigma of gauss kernel with which grid should be smoothed prior to re-distancing.
+ */
+ void run_gauss_filter(const size_t sigma)
+ {
+ size_t n = count_gauss_repitition(sigma);
+ for (size_t i = 0; i < n; i++)
+ {
+ gauss_filter<Phi_0_temp, Phi_nplus1_temp>(g_temp);
+ copy_gridTogrid<Phi_nplus1_temp, Phi_0_temp>(g_temp, g_temp);
+ }
+ }
+
+ /** @brief Run one timestep of re-distancing and compute Phi_n+1.
+ *
+ * @param phi_n Phi value on current node and current time.
+ * @param phi_n_magnOfGrad Gradient magnitude of current Phi from upwinding FD.
+ * @param dt Time step.
+ * @param sign_phi_n Sign of the current Phi, should be the smooth sign.
+ *
+ * @return Phi_n+1 which is the Phi of the next time step on current node.
+ *
+ */
+ double get_phi_nplus1(double phi_n, double phi_n_magnOfGrad, double dt, double sign_phi_n)
+ {
+ double step = dt * sign_phi_n * (1 - phi_n_magnOfGrad); // <- original Sussman
+ if (step > 10)
+ {
+ std::cout << "phi_n_magnOfGrad = " << phi_n_magnOfGrad << ", step = " << step
+ << ", skip to prevent exploding peaks." << std::endl;
+ step = 0;
+ }
+ return phi_n + step;
+ }
+
+ /** @brief Go one re-distancing time-step on the whole grid.
+ *
+ * @param grid Internal temporary grid.
+ */
+ void go_one_redistancing_step_whole_grid(g_temp_type &grid)
+ {
+ grid.template ghost_get<Phi_0_temp, Phi_nplus1_temp, Phi_grad_temp, Phi_magnOfGrad_temp>();
+ double spacing_x = grid.getSpacing()[0];
+ auto dom = grid.getDomainIterator();
+ while (dom.isNext())
+ {
+ auto key = dom.get();
+ const double phi_n = grid.template get<Phi_0_temp>(key);
+ const double phi_n_magnOfGrad = grid.template get<Phi_magnOfGrad_temp>(key);
+ double epsilon = phi_n_magnOfGrad * spacing_x;
+ grid.template get<Phi_nplus1_temp>(key) = get_phi_nplus1(phi_n, phi_n_magnOfGrad, time_step,
+ smooth_S(phi_n, epsilon));
+ ++dom;
+ }
+ }
+
+ /** @brief Updates Phi_n with the new Phi_n+1 and recomputes the gradients.
+ *
+ * @param grid Internal temporary grid.
+ */
+ void update_grid(g_temp_type &grid)
+ {
+ copy_gridTogrid<Phi_nplus1_temp, Phi_0_temp>(grid, grid); // Update Phi_0
+ get_upwind_gradient<Phi_0_temp, Phi_0_sign_temp, Phi_grad_temp>(grid);
+ get_gradient_magnitude<Phi_grad_temp, Phi_magnOfGrad_temp>(grid);
+ }
+
+ /** @brief Checks if a node lays within the narrow band around the interface.
+ *
+ * @param Phi Value of Phi at that specific node.
+ *
+ * @return True, if node lays within nb., false, if the distance to the interface is > kappa.
+ */
+ bool lays_inside_NB(double Phi)
+ {
+ return (abs(Phi) <= kappa);
+ }
+
+ /** @brief Checks how far current solution is from fulfilling the user-defined convergence criteria.
+ *
+ * @param grid Internal temporary grid.
+ *
+ * @return Total residual (1 - phi_gradient_magnitude) and total change from the last time step,
+ * both normalized by the number of grid nodes in the narrow band.
+ */
+ DistFromSol get_residual_and_change_NB(g_temp_type &grid)
+ {
+ double total_residual = 0;
+ double total_change = 0;
+ double total_points_in_nb = 0;
+ auto dom = grid.getDomainIterator();
+ while (dom.isNext())
+ {
+ auto key = dom.get();
+ if (lays_inside_NB(grid.template get<Phi_nplus1_temp>(key)))
+ {
+ total_points_in_nb += 1.0;
+ total_residual += abs(grid.template get<Phi_magnOfGrad_temp>(key) - 1);
+ total_change += abs(grid.template get<Phi_nplus1_temp>(key) - grid.template get<Phi_0_temp>(key));
+ }
+ ++dom;
+ }
+ auto &v_cl = create_vcluster();
+ v_cl.sum(total_points_in_nb);
+ v_cl.sum(total_residual);
+ v_cl.sum(total_change);
+ v_cl.execute();
+ return {total_residual / total_points_in_nb, total_change / total_points_in_nb};
+ }
+
+ /** @brief Prints out the iteration number, residual and change of the current re-distancing iteration
+ *
+ * @param grid Internal temporary grid.
+ * @param iter Current re-distancing iteration.
+ *
+ * @return Total residual (1 - phi_gradient_magnitude) and total change from the last time step,
+ * both normalized by the number of grid nodes in the narrow band.
+ */
+ void print_out_iteration_change_residual(g_temp_type &grid, size_t iter)
+ {
+ DistFromSol distFromSol = get_residual_and_change_NB(grid);
+ auto &v_cl = create_vcluster();
+ if (v_cl.rank() == 0)
+ {
+ if (iter == 0)
+ {
+ std::cout << "Iteration,Change,Residual" << std::endl;
+ }
+ std::cout << iter << "," << distFromSol.change << "," << distFromSol.residual << std::endl;
+ }
+ }
+
+ /** @brief Checks steady-state is reached in the narrow band.
+ *
+ * @details Checks if change and/or residual between 2 iterations smaller than user-defined convolution tolerance.
+ *
+ * @param grid Internal temporary grid.
+ *
+ * @return True, if steady-state reached, else false.
+ *
+ * @see Conv_tol_change, Conv_tol_residual
+ */
+ bool steady_state_NB(g_temp_type &grid)
+ {
+ bool steady_state = false;
+ DistFromSol distFromSol = get_residual_and_change_NB(grid);
+ if (redistOptions.convTolChange.check && redistOptions.convTolResidual.check)
+ {
+ steady_state = (distFromSol.change <= redistOptions.convTolChange.value &&
+ distFromSol.residual <= redistOptions.convTolResidual.value);
+ }
+ else
+ {
+ if (redistOptions.convTolChange.check)
+ {
+ steady_state = (distFromSol.change <= redistOptions.convTolChange.value);
+ } // Use the normalized total change between two iterations in the narrow bands steady-state criterion
+ if (redistOptions.convTolResidual.check)
+ {
+ steady_state = (distFromSol.residual <= redistOptions.convTolResidual.value);
+ } // Use the normalized total residual of phi compared to SDF in the narrow bands steady-state criterion
+ }
+ return steady_state;
+ }
+
+ /** @brief Runs Sussman re-distancing on the internal temporary grid.
+ *
+ * @details The number of iterations is minimum redistOptions.min_iter iterations and finishes when either
+ * steady-state or redist_options.max_iter is reached. The steady-state convergence accuracy depends on the
+ * user defined redist_options.convTolChange.value and redist_options.convTolResidual.value, respectively.
+ *
+ * @param grid Internal temporary grid.
+ *
+ */
+ void iterative_redistancing(g_temp_type &grid)
+ {
+ for (size_t i = 0; i <= redistOptions.max_iter; i++)
+ {
+ go_one_redistancing_step_whole_grid(grid);
+
+ if (redistOptions.print_current_iterChangeResidual)
+ {
+ print_out_iteration_change_residual(grid, i);
+ }
+
+ if (i >= redistOptions.min_iter)
+ {
+ if (i % redistOptions.interval_check_convergence ==
+ 0) // after some iterations check if steady state is reached in the narrow band
+ {
+ if (steady_state_NB(grid))
+ {
+ if (redistOptions.print_steadyState_iter)
+ {
+ auto &v_cl = create_vcluster();
+ if (v_cl.rank() == 0)
+ {
+ std::cout << "Steady state reached at iteration: " << i << std::endl;
+ std::cout << "Final_Iteration,Change,Residual" << std::endl;
+ }
+ print_out_iteration_change_residual(grid, i);
+ }
+ update_grid(grid); // Update Phi
+ break;
+ }
+ auto &v_cl = create_vcluster();
+ }
+ }
+ update_grid(grid);
+ }
+ }
+};
+
+#endif //REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP