Commit e5c82bd7 authored by jstark's avatar jstark
Browse files

Generalizing types in analytical signed distance functions. Adding a fast unit...

Generalizing types in analytical signed distance functions. Adding a fast unit test for Sussman redistancing.
parent d0e2c5cd
......@@ -29,8 +29,14 @@
template <typename point_type, typename radius_type>
bool inside_sphere(point_type coords, radius_type radius, double center_x=0, double center_y=0, double center_z=0)
{
const double EPSILON = std::numeric_limits<double>::epsilon();
const double X = coords.get(0), Y = coords.get(1), Z = coords.get(2);
typedef typename std::remove_const_t<std::remove_reference_t<decltype(coords.get(0))>> space_type;
if(!(std::is_same<space_type, radius_type>::value))
{
std::cout << "Radius type and space type of grid must be the same! Aborting..." << std::endl;
abort();
}
const space_type EPSILON = std::numeric_limits<space_type>::epsilon();
const space_type X = coords.get(0), Y = coords.get(1), Z = coords.get(2);
return (X - center_x) * (X - center_x)
+ (Y - center_y) * (Y - center_y)
+ (Z - center_z) * (Z - center_z)
......@@ -57,13 +63,18 @@ bool inside_sphere(point_type coords, radius_type radius, double center_x=0, dou
template <size_t Phi_0, typename grid_type, typename radius_type>
void init_grid_with_sphere(grid_type & grid, radius_type radius, double center_x=0, double center_y=0, double center_z=0)
{
if(!(std::is_same<typename grid_type::stype, radius_type>::value))
{
std::cout << "Radius type and space type of grid must be the same! Aborting..." << std::endl;
abort();
}
// 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();
Point<grid_type::dims, double> coords = grid.getPos(key);
Point<grid_type::dims, typename grid_type::stype> coords = grid.getPos(key);
if (inside_sphere(coords, radius, center_x, center_y, center_z))
{
......
......@@ -73,8 +73,8 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
// get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm forward FD");
// write_lnorms_to_file(N, lNorms, "l_norms_FDfwd", "./");
std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
......@@ -134,8 +134,8 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
// get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm backward FD");
// write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
......@@ -195,8 +195,8 @@ BOOST_AUTO_TEST_SUITE(FDOrder1TestSuite)
// get_absolute_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
get_relative_error<NumericalGradient, AnalyticalGradient, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
BOOST_CHECK_MESSAGE(lNorms.l2 < l2_norms[count] + 0.00001 + EPSILON, "Checking L2-norm central FD");
// write_lnorms_to_file(N, lNorms, "l_norms_FDbwd", "./");
std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
......
......@@ -74,11 +74,11 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
// std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
write_lnorms_to_file(N, lNorms, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
lNorms.write_to_file(N, 6, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
"./");
// if(N==128) g_dist.write("grid_gaussian_upwindGradient_N" + std::to_string(N), FORMAT_BINARY);
......@@ -145,11 +145,11 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
// std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
// write_lnorms_to_file(N, lNorms, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
// lNorms.write_to_file(N, 6, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
// "./");
// if(N==128) g_dist.write("grid_gaussian_upwindGradient_N" + std::to_string(N), FORMAT_BINARY);
......@@ -213,11 +213,11 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
// std::cout << N << ", " << lNorms.l2 << ", " << lNorms.linf << std::endl;
// write_lnorms_to_file(N, lNorms, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
// lNorms.write_to_file(N, 6, "l_norms_upwindGrad_convOrder_" + std::to_string(convergence_order), ""
// "./");
// if(N==128) g_dist.write("grid_gaussian_upwindGradient_N" + std::to_string(N), FORMAT_BINARY);
......@@ -280,9 +280,10 @@ BOOST_AUTO_TEST_SUITE(UpwindGradientTestSuite)
// Get the error between analytical and numerical solution
get_relative_error<df_upwind, df_gaussian, Error>(g_dist);
L_norms lNorms;
lNorms = get_l_norms_grid<Error>(g_dist);
// write_lnorms_to_file(N, lNorms, "l_norms_upwindGrad_3D_convOrder_" + std::to_string(convergence_order), ""
LNorms<double> lNorms;
lNorms.get_l_norms_grid<Error>(g_dist);
// lNorms.write_to_file(N, 6, "l_norms_upwindGrad_3D_convOrder_" + std::to_string(convergence_order)
// , ""
// "./");
// g_dist.write("grid_gaussian_upwindGradient_N" + std::to_string(N) + "_order" + std::to_string
// (convergence_order), FORMAT_BINARY);
......
......@@ -39,11 +39,17 @@
* @return Double variable that contains the exact solution for the signed distance function of a given point in a
* sphere of given radius, where the SDF has positive values inside and negative values outside the sphere.
*/
template <typename point_type, typename radius_type, typename center_type>
double get_analytic_sdf_sphere(point_type coords, radius_type radius,
center_type center_x=0, center_type center_y=0, center_type center_z=0)
template <typename point_type, typename space_type>
space_type get_analytic_sdf_sphere(point_type coords, space_type radius,
space_type center_x=0, space_type center_y=0, space_type center_z=0)
{
const double X = coords.get(0), Y = coords.get(1), Z = coords.get(2);
typedef typename std::remove_const_t<std::remove_reference_t<decltype(coords.get(0))>> coord_type;
if(!(std::is_same<space_type, coord_type>::value))
{
std::cout << "Radius-, Center- and Space-type of grid must be the same! Aborting..." << std::endl;
abort();
}
const space_type X = coords.get(0), Y = coords.get(1), Z = coords.get(2);
return (radius -
sqrt((X - center_x) * (X - center_x)
+ (Y - center_y) * (Y - center_y)
......@@ -64,15 +70,21 @@ double get_analytic_sdf_sphere(point_type coords, radius_type radius,
* @param center_z Double z-coordinate of sphere center.
*/
template <size_t SDF_exact, typename grid_type, typename radius_type, typename center_type>
void init_analytic_sdf_sphere(grid_type & grid, radius_type radius, center_type center_x=0, center_type center_y=0,
center_type center_z=0)
template <size_t SDF_exact, typename grid_type, typename space_type>
void init_analytic_sdf_sphere(grid_type & grid, space_type radius, space_type center_x=0, space_type center_y=0,
space_type center_z=0)
{
if(!(std::is_same<typename grid_type::stype, space_type>::value))
{
std::cout << "Radius-, Center- and Space-type of grid must be the same! Aborting..." << std::endl;
abort();
}
auto dom = grid.getDomainIterator();
while(dom.isNext())
{
auto key = dom.get();
Point<grid_type::dims, double> coords = grid.getPos(key);
Point<grid_type::dims, typename grid_type::stype> coords = grid.getPos(key);
grid.template getProp<SDF_exact>(key) = get_analytic_sdf_sphere(coords, radius, center_x,
center_y, center_z);
++dom;
......@@ -97,11 +109,17 @@ void init_analytic_sdf_sphere(grid_type & grid, radius_type radius, center_type
* @return Double variable that contains the exact solution for the signed distance function of a given point in a
* sphere of given radius, where the SDF has positive values inside and negative values outside the sphere.
*/
template <typename point_type, typename radius_type, typename center_type>
double get_analytic_sdf_circle(point_type coords, radius_type radius,
center_type center_x=0, center_type center_y=0)
template <typename point_type, typename space_type>
space_type get_analytic_sdf_circle(point_type coords, space_type radius,
space_type center_x=0, space_type center_y=0)
{
const double X = coords.get(0), Y = coords.get(1);
typedef typename std::remove_const_t<std::remove_reference_t<decltype(coords.get(0))>> coord_type;
if(!(std::is_same<space_type, coord_type>::value))
{
std::cout << "Radius-, Center- and Space-type of grid must be the same! Aborting..." << std::endl;
abort();
}
const space_type X = coords.get(0), Y = coords.get(1);
return (radius -
sqrt((X - center_x) * (X - center_x)
+ (Y - center_y) * (Y - center_y)));
......@@ -120,14 +138,19 @@ double get_analytic_sdf_circle(point_type coords, radius_type radius,
* @param center_x X-coordinate of the circle center.
* @param center_y Y-coordinate of the circle center.
*/
template <size_t SDF_exact, typename grid_type, typename radius_type, typename center_type>
void init_analytic_sdf_circle(grid_type & grid, radius_type radius, center_type center_x=0, center_type center_y=0)
template <size_t SDF_exact, typename grid_type, typename space_type>
void init_analytic_sdf_circle(grid_type & grid, space_type radius, space_type center_x=0, space_type center_y=0)
{
if(!(std::is_same<typename grid_type::stype, space_type>::value))
{
std::cout << "Radius-, Center- and Space-type of grid must be the same! Aborting..." << std::endl;
abort();
}
auto dom = grid.getDomainIterator();
while(dom.isNext())
{
auto key = dom.get();
Point<grid_type::dims, double> coords = grid.getPos(key);
Point<grid_type::dims, typename grid_type::stype> coords = grid.getPos(key);
grid.template getProp<SDF_exact>(key) = get_analytic_sdf_circle(coords, radius, center_x,
center_y);
++dom;
......
......@@ -122,9 +122,9 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
// vd_narrow_band.write("vd_error_N" + std::to_string(N), FORMAT_BINARY);
// vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
// Compute the L_2- and L_infinity-norm and save to file
L_norms lNorms_vd;
lNorms_vd = get_l_norms_vector<Error_vd>(vd_narrow_band);
write_lnorms_to_file(N, lNorms_vd, "l_norms_", "./");
LNorms<phi_type> lNorms_vd;
lNorms_vd.get_l_norms_vector<Error_vd>(vd_narrow_band);
lNorms_vd.write_to_file(N, 6, "l_norms_", "./");
}
}
......@@ -209,9 +209,9 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
// std::to_string(order), FORMAT_BINARY);
// vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
// Compute the L_2- and L_infinity-norm and save to file
L_norms lNorms_vd;
lNorms_vd = get_l_norms_vector<Error_vd>(vd_narrow_band);
write_lnorms_to_file(N, lNorms_vd, "l_norms_vd_absError_8p_order" + std::to_string(order), "./");
LNorms<phi_type> lNorms_vd;
lNorms_vd.get_l_norms_vector<Error_vd>(vd_narrow_band);
lNorms_vd.write_to_file(N, 6, "l_norms_vd_absError_8p_order" + std::to_string(order), "./");
// switch(order)
// {
......
......@@ -27,13 +27,6 @@
#include "level_set/redistancing_Sussman/HelpFunctions.hpp" // for printing to_string_with_precision
/**@brief Structure that bundles the two variables for L_2 and L_infinity norm.
*/
struct L_norms
{
double l2; ///< Double variable that is supposed to contain the L_2 norm.
double linf; ///< Double variable that is supposed to contain the L_infinity norm.
};
/**@brief At each grid node, the absolute error is computed and stored in another property.
*
* @details The absolute error is computed as:
......@@ -50,7 +43,16 @@ template <size_t PropNumeric, size_t PropAnalytic, size_t Error, typename gridty
void get_absolute_error(gridtype & grid)
{
auto dom = grid.getDomainIterator();
while(dom.isNext())
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
grid.template get<PropNumeric>(dom.get()))>> numeric_type;
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
grid.template get<PropAnalytic>(dom.get()))>> analytic_type;
if(!(std::is_same<numeric_type, analytic_type>::value))
{
std::cout << "Type of numerical and analytical solution must be the same! Aborting..." << std::endl;
abort();
} while(dom.isNext())
{
auto key = dom.get();
grid.template getProp<Error> (key) = abs(grid.template getProp<PropAnalytic> (key) - (grid.template getProp<PropNumeric> (key)));
......@@ -72,8 +74,18 @@ void get_absolute_error(gridtype & grid)
template <size_t PropNumeric, size_t PropAnalytic, size_t Error, typename gridtype>
void get_relative_error(gridtype & grid)
{
const double EPSILON = std::numeric_limits<double>::epsilon();
auto dom = grid.getDomainIterator();
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
grid.template get<PropNumeric>(dom.get()))>> numeric_type;
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
grid.template get<PropAnalytic>(dom.get()))>> analytic_type;
if(!(std::is_same<numeric_type, analytic_type>::value))
{
std::cout << "Type of numerical and analytical solution must be the same! Aborting..." << std::endl;
abort();
}
while(dom.isNext())
{
auto key = dom.get();
......@@ -82,6 +94,21 @@ void get_relative_error(gridtype & grid)
++dom;
}
}
/**@brief Class for computing relative/absolute errors and l2/l_infinity norm for distributed grids and vectors
*
* @tparam lnorm_type Desired return type for l-norm.
*/
template <typename lnorm_type>
class LNorms
{
public:
LNorms(){};
// Member variables
lnorm_type l2;
lnorm_type linf;
// Member functions
/**@brief Computes the L_2 and L_infinity norm on the basis of the precomputed error on a grid.
*
* @tparam Error Index of grid property that contains the error.
......@@ -89,27 +116,30 @@ void get_relative_error(gridtype & grid)
* @param grid Input OpenFPM grid. Can be of any dimension.
* @return Object of type L_norms that contains #L_norms::l2 and #L_norms::linf.
*/
template <size_t Error, typename gridtype>
L_norms get_l_norms_grid(gridtype & grid)
{
double maxError = 0;
double sumErrorSq = 0;
auto dom = grid.getDomainIterator();
while(dom.isNext())
template <size_t Error, typename gridtype>
void get_l_norms_grid(gridtype & grid)
{
auto key = dom.get();
sumErrorSq += grid.template getProp<Error> (key) * grid.template getProp<Error> (key);
if (grid.template getProp<Error> (key) > maxError) maxError = grid.template getProp<Error> (key); // update maxError
++dom;
auto dom = grid.getDomainIterator();
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
grid.template get<Error>(dom.get()))>> error_type;
error_type maxError = 0;
error_type sumErrorSq = 0;
while(dom.isNext())
{
auto key = dom.get();
sumErrorSq += grid.template getProp<Error> (key) * grid.template getProp<Error> (key);
if (grid.template getProp<Error> (key) > maxError) maxError = grid.template getProp<Error> (key); // update maxError
++dom;
}
auto &v_cl = create_vcluster();
v_cl.sum(sumErrorSq);
v_cl.max(maxError);
v_cl.execute();
l2 = (lnorm_type) sqrt( sumErrorSq / (error_type)grid.size());
linf = (lnorm_type) maxError;
}
auto &v_cl = create_vcluster();
v_cl.sum(sumErrorSq);
v_cl.max(maxError);
v_cl.execute();
double l2 = sqrt( sumErrorSq / (double)grid.size());
double linf = maxError;
return {l2, linf};
}
/**@brief Computes the L_2 and L_infinity norm on the basis of the precomputed error on a particle vector.
*
* @tparam Error Index of grid property that contains the error.
......@@ -117,30 +147,32 @@ L_norms get_l_norms_grid(gridtype & grid)
* @param vd Input particle vector.
* @return Object of type L_norms that contains #L_norms::l2 and #L_norms::linf.
*/
template <size_t Error, typename vectortype>
L_norms get_l_norms_vector(vectortype & vd)
{
double maxError = 0;
double sumErrorSq = 0;
auto dom = vd.getDomainIterator();
double count = 0;
while(dom.isNext())
template <size_t Error, typename vectortype>
void get_l_norms_vector(vectortype & vd)
{
auto key = dom.get();
sumErrorSq += vd.template getProp<Error> (key) * vd.template getProp<Error> (key);
if (vd.template getProp<Error> (key) > maxError) maxError = vd.template getProp<Error> (key); // update maxError
++dom;
++count;
auto dom = vd.getDomainIterator();
typedef typename std::remove_const_t<std::remove_reference_t<decltype(
vd.template getProp<Error>(dom.get()))>> error_type;
error_type maxError = 0;
error_type sumErrorSq = 0;
int count = 0;
while(dom.isNext())
{
auto key = dom.get();
sumErrorSq += vd.template getProp<Error> (key) * vd.template getProp<Error> (key);
if (vd.template getProp<Error> (key) > maxError) maxError = vd.template getProp<Error> (key); // update maxError
++dom;
++count;
}
auto &v_cl = create_vcluster();
v_cl.sum(sumErrorSq);
v_cl.sum(count);
v_cl.max(maxError);
v_cl.execute();
l2 = (lnorm_type) sqrt( sumErrorSq / (error_type)count);
linf = (lnorm_type) maxError;
}
auto &v_cl = create_vcluster();
v_cl.sum(sumErrorSq);
v_cl.sum(count);
v_cl.max(maxError);
v_cl.execute();
double l2 = sqrt( sumErrorSq / (double)count);
double linf = maxError;
return {l2, linf};
}
/**@brief Writes the N (number of grid points on a square grid) and L-norms as strings to a csv-file.
*
* @param N Size_t variable that contains the grid size in number of grid points in one dimension for an NxN(xN) grid
......@@ -149,22 +181,33 @@ L_norms get_l_norms_vector(vectortype & vd)
* written to.
* @param path_output Std::string containing the path where the output csv file should be saved.
*/
static void write_lnorms_to_file(size_t N, L_norms l_norms, std::string filename, std::string path_output)
{
auto &v_cl = create_vcluster();
if (v_cl.rank() == 0)
void write_to_file(const size_t N, const int precision,
const std::string & filename, const std::string & path_output)
{
std::string path_output_lnorm = path_output + "/" + filename + ".csv";
create_file_if_not_exist(path_output_lnorm);
std::ofstream l_out;
l_out.open(path_output_lnorm, std::ios_base::app); // append instead of overwrite
l_out << std::to_string(N)
<< ',' << to_string_with_precision(l_norms.l2, 15)
<< ',' << to_string_with_precision(l_norms.linf, 15) << std::endl;
l_out.close();
auto &v_cl = create_vcluster();
if (v_cl.rank() == 0)
{
std::string path_output_lnorm = path_output + "/" + filename + ".csv";
create_file_if_not_exist(path_output_lnorm);
std::ofstream l_out;
l_out.open(path_output_lnorm, std::ios_base::app); // append instead of overwrite
l_out << std::to_string(N)
<< ',' << to_string_with_precision(l2, precision)
<< ',' << to_string_with_precision(linf, precision) << std::endl;
l_out.close();
}
}
}
private:
};
......
//
// Created by jstark on 04.01.22.
//
#define BOOST_TEST_DYN_LINK
//#define BOOST_TEST_MAIN // in only one cpp file
#include <boost/test/unit_test.hpp>
// Include redistancing files
#include "util/PathsAndFiles.hpp"
#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
#include "level_set/redistancing_Sussman/NarrowBand.hpp"
// Include header files for testing
#include "Draw/DrawSphere.hpp"
#include "l_norms/LNorms.hpp"
#include "analytical_SDF/AnalyticalSDF.hpp"
BOOST_AUTO_TEST_SUITE(RedistancingSussmanFastTestSuite)
BOOST_AUTO_TEST_CASE(RedistancingSussman_unit_sphere_fast_double_test)
{
typedef double phi_type;
typedef double space_type;
const phi_type EPSILON = std::numeric_limits<phi_type>::epsilon();
const size_t grid_dim = 3;
// some indices
const size_t x = 0;
const size_t y = 1;
const size_t z = 2;
const size_t Phi_0_grid = 0;
const size_t SDF_sussman_grid = 1;
const size_t SDF_exact_grid = 2;
const size_t Error_grid = 3;
size_t N = 32;
const size_t sz[grid_dim] = {N, N, N};
const space_type radius = 1.0;
const space_type box_lower = -2.0;
const space_type box_upper = 2.0;
Box<grid_dim, space_type> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
Ghost<grid_dim, long int> ghost(0);
typedef aggregate<phi_type, phi_type, phi_type, phi_type> props;
typedef grid_dist_id<grid_dim, space_type, props > grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
const space_type center[grid_dim] = {(box_upper+box_lower)/(space_type)2,
(box_upper+box_lower)/(space_type)2,
(box_upper+box_lower)/(space_type)2};
init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
Redist_options<phi_type> redist_options;
redist_options.min_iter = 1e3;
redist_options.max_iter = 1e3;
redist_options.convTolChange.check = false;
redist_options.convTolResidual.check = false;
redist_options.interval_check_convergence = 1e3;
redist_options.width_NB_in_grid_points = 8;
redist_options.print_current_iterChangeResidual = true;
redist_options.print_steadyState_iter = true;
RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options); // Instantiation of
std::cout << "Automatically found timestep is " << redist_obj.get_time_step() << std::endl;
// Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
redist_obj.run_redistancing<Phi_0_grid, SDF_sussman_grid>();
// Compute exact signed distance function at each grid point
init_analytic_sdf_sphere<SDF_exact_grid>(g_dist, radius, center[x], center[y], center[z]);
// Compute the absolute error between analytical and numerical solution at each grid point
get_absolute_error<SDF_sussman_grid, SDF_exact_grid, Error_grid>(g_dist);
g_dist.write("/Users/jstark/Desktop/unit_testing/grid_post_redistancing_fast_test", FORMAT_BINARY);
/////////////////////////////////////////////////////////////////////////////////////////////
// Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
// interface)
size_t bc[grid_dim] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
typedef aggregate<phi_type> props_nb;
typedef vector_dist<grid_dim, space_type, props_nb> vd_type;
Ghost<grid_dim, space_type> ghost_vd(0);
vd_type vd_narrow_band(0, box, bc, ghost_vd);
vd_narrow_band.setPropNames({"error"});
size_t narrow_band_width = 8;
NarrowBand<grid_in_type, phi_type> narrowBand(g_dist, narrow_band_width); // Instantiation of NarrowBand class
const size_t Error_vd = 0;
narrowBand.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
vd_narrow_band);
vd_narrow_band.write("/Users/jstark/Desktop/unit_testing/vd_narrow_band_fast_test", FORMAT_BINARY);
// Compute the L_2- and L_infinity-norm and save to file
LNorms<phi_type> lNorms_vd;