Adding redistancing Sussmann

src/level_set/redistancing_Sussman/ComputeGradient.hpp

+//
+// 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

src/level_set/redistancing_Sussman/GaussFilter.hpp

+//
+// 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

src/level_set/redistancing_Sussman/HelpFunctions.hpp

+//
+// 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 &phi; 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

src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp

+//
+// 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

src/level_set/redistancing_Sussman/NarrowBand.hpp

+//
+// 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

src/level_set/redistancing_Sussman/RedistancingSussman.hpp

+//
+// 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 &phi; 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 &phi; 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