diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8c07f8e59c0a175abf5f02cbd536c92937a3b995..da44646be62bf12f260248645a60d510b57dfd1d 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -33,74 +33,74 @@ if ( CUDA_ON_BACKEND STREQUAL "HIP" AND HIP_FOUND )
hip_add_executable(numerics ${OPENFPM_INIT_FILE} ${CUDA_SOURCES}
- OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
- FiniteDifference/FD_Solver_test.cpp
- FiniteDifference/FD_op_Tests.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test_temporal.cpp
- DCPSE/tests/Dcpse_unit_tests.cpp
- DCPSE/tests/DcpseRhs_unit_tests.cpp
- DCPSE/tests/MonomialBasis_unit_tests.cpp
- DCPSE/tests/Support_unit_tests.cpp
- DCPSE/tests/Vandermonde_unit_tests.cpp
+# OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
+# FiniteDifference/FD_Solver_test.cpp
+# FiniteDifference/FD_op_Tests.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test_temporal.cpp
+# DCPSE/tests/Dcpse_unit_tests.cpp
+# DCPSE/tests/DcpseRhs_unit_tests.cpp
+# DCPSE/tests/MonomialBasis_unit_tests.cpp
+# DCPSE/tests/Support_unit_tests.cpp
+# DCPSE/tests/Vandermonde_unit_tests.cpp
main.cpp
- Matrix/SparseMatrix_unit_tests.cpp
- interpolation/interpolation_unit_tests.cpp
- Vector/Vector_unit_tests.cpp
- Solvers/petsc_solver_unit_tests.cpp
- FiniteDifference/FDScheme_unit_tests.cpp
- FiniteDifference/eq_unit_test_3d.cpp
- FiniteDifference/eq_unit_test.cpp
- FiniteDifference/tests/Eno_Weno_unit_test.cpp
- FiniteDifference/tests/Upwind_gradient_unit_test.cpp
- FiniteDifference/tests/FD_simple_unit_test.cpp
- Operators/Vector/vector_dist_operators_unit_tests.cpp
- Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
- ../../src/lib/pdata.cpp
- BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
+# Matrix/SparseMatrix_unit_tests.cpp
+# interpolation/interpolation_unit_tests.cpp
+# Vector/Vector_unit_tests.cpp
+# Solvers/petsc_solver_unit_tests.cpp
+# FiniteDifference/FDScheme_unit_tests.cpp
+# FiniteDifference/eq_unit_test_3d.cpp
+# FiniteDifference/eq_unit_test.cpp
+# FiniteDifference/tests/Eno_Weno_unit_test.cpp
+# FiniteDifference/tests/Upwind_gradient_unit_test.cpp
+# FiniteDifference/tests/FD_simple_unit_test.cpp
+# Operators/Vector/vector_dist_operators_unit_tests.cpp
+# Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
+# ../../src/lib/pdata.cpp
+# BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
# level_set/closest_point/closest_point_unit_tests.cpp
-# level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
-# level_set/redistancing_Sussman/tests/convergence_test.cpp
+ level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
+ level_set/redistancing_Sussman/tests/convergence_test.cpp
)
else()
add_executable(numerics ${OPENFPM_INIT_FILE} ${CUDA_SOURCES}
- OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
- FiniteDifference/FD_Solver_test.cpp
- FiniteDifference/FD_op_Tests.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cpp
- DCPSE/DCPSE_op/tests/DCPSE_op_test_temporal.cpp
- DCPSE/tests/Dcpse_unit_tests.cpp
- DCPSE/tests/DcpseRhs_unit_tests.cpp
- DCPSE/tests/MonomialBasis_unit_tests.cpp
- DCPSE/tests/Support_unit_tests.cpp
- DCPSE/tests/Vandermonde_unit_tests.cpp
+# OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
+# FiniteDifference/FD_Solver_test.cpp
+# FiniteDifference/FD_op_Tests.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cpp
+# DCPSE/DCPSE_op/tests/DCPSE_op_test_temporal.cpp
+# DCPSE/tests/Dcpse_unit_tests.cpp
+# DCPSE/tests/DcpseRhs_unit_tests.cpp
+# DCPSE/tests/MonomialBasis_unit_tests.cpp
+# DCPSE/tests/Support_unit_tests.cpp
+# DCPSE/tests/Vandermonde_unit_tests.cpp
main.cpp
- Matrix/SparseMatrix_unit_tests.cpp
- interpolation/interpolation_unit_tests.cpp
- Vector/Vector_unit_tests.cpp
- Solvers/petsc_solver_unit_tests.cpp
- FiniteDifference/FDScheme_unit_tests.cpp
- FiniteDifference/eq_unit_test_3d.cpp
- FiniteDifference/eq_unit_test.cpp
- FiniteDifference/tests/Eno_Weno_unit_test.cpp
- FiniteDifference/tests/Upwind_gradient_unit_test.cpp
- FiniteDifference/tests/FD_simple_unit_test.cpp
- Operators/Vector/vector_dist_operators_unit_tests.cpp
- Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
- ../../src/lib/pdata.cpp
- BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
+# Matrix/SparseMatrix_unit_tests.cpp
+# interpolation/interpolation_unit_tests.cpp
+# Vector/Vector_unit_tests.cpp
+# Solvers/petsc_solver_unit_tests.cpp
+# FiniteDifference/FDScheme_unit_tests.cpp
+# FiniteDifference/eq_unit_test_3d.cpp
+# FiniteDifference/eq_unit_test.cpp
+# FiniteDifference/tests/Eno_Weno_unit_test.cpp
+# FiniteDifference/tests/Upwind_gradient_unit_test.cpp
+# FiniteDifference/tests/FD_simple_unit_test.cpp
+# Operators/Vector/vector_dist_operators_unit_tests.cpp
+# Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
+# ../../src/lib/pdata.cpp
+# BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp
# level_set/closest_point/closest_point_unit_tests.cpp
-# level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
-# level_set/redistancing_Sussman/tests/convergence_test.cpp
+ level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
+ level_set/redistancing_Sussman/tests/convergence_test.cpp
)
set_property(TARGET numerics PROPERTY CUDA_ARCHITECTURES OFF)
diff --git a/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
index 2678679613a40d2374cdd4c295da5a9fe23b3a64..ee47be5fd6ee46a3718c349b1723c076e5e28152 100644
--- a/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
+++ b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
@@ -26,9 +26,9 @@
* @return Time step.
*/
template <typename grid_type>
-double get_time_step_CFL(grid_type &grid)
+typename grid_type::stype get_time_step_CFL(grid_type &grid)
{
- double sum = 0;
+ typename grid_type::stype sum = 0.0;
for (size_t d = 0; d < grid_type::dims; d++)
{
sum += 1.0 / (grid.spacing(d) * grid.spacing(d));
@@ -47,9 +47,9 @@ double get_time_step_CFL(grid_type &grid)
* @return Time step.
*/
template <typename grid_type>
-double get_time_step_CFL(grid_type & grid, double u [grid_type::dims], double C)
+typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type::stype u [grid_type::dims], typename grid_type::stype C)
{
- double sum = 0;
+ typename grid_type::stype sum = 0;
for (size_t d = 0; d < grid_type::dims; d++)
{
sum += u[d] / grid.spacing(d);
@@ -66,9 +66,9 @@ double get_time_step_CFL(grid_type & grid, double u [grid_type::dims], double C)
* @return Time step.
*/
template <typename grid_type>
-double get_time_step_CFL(grid_type & grid, double u, double C)
+typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type::stype u, typename grid_type::stype C)
{
- double sum = 0;
+ typename grid_type::stype sum = 0;
for (size_t d = 0; d < grid_type::dims; d++)
{
sum += u / grid.spacing(d);
@@ -149,9 +149,9 @@ void copy_gridTogrid(const grid_source_type & grid_sc, grid_dest_type & grid_ds,
* @return Double variable that contains the value of the biggest spacing.
*/
template <typename grid_type>
-double get_biggest_spacing(grid_type & grid)
+typename grid_type::stype get_biggest_spacing(grid_type & grid)
{
- double h_max = 0;
+ typename grid_type::stype h_max = 0;
for (size_t d = 0; d < grid_type::dims; d++)
{
if (grid.spacing(d) > h_max) h_max = grid.spacing(d);
@@ -167,9 +167,9 @@ double get_biggest_spacing(grid_type & grid)
* @return Double variable that contains the value of the smallest spacing.
*/
template <typename grid_type>
-double get_smallest_spacing(grid_type & grid)
+typename grid_type::stype get_smallest_spacing(grid_type & grid)
{
- double spacing [grid_type::dims];
+ typename grid_type::stype spacing [grid_type::dims];
for (size_t d = 0; d < grid_type::dims; d++)
{
spacing[d] = grid.spacing(d);
@@ -188,10 +188,10 @@ double get_smallest_spacing(grid_type & 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)
+template <size_t Prop1, size_t Prop2, typename prop_type, typename grid_type>
+prop_type average_difference(grid_type & grid)
{
- double total_diff = 0;
+ prop_type total_diff = 0;
auto dom = grid.getDomainIterator();
while (dom.isNext())
{
@@ -209,10 +209,10 @@ double average_difference(grid_type & 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)
+template <size_t Prop, typename prop_type, typename grid_type>
+prop_type get_max_val(grid_type & grid)
{
- double max_value = std::numeric_limits<double>::lowest();
+ prop_type max_value = std::numeric_limits<prop_type>::lowest();
auto dom = grid.getDomainIterator();
while(dom.isNext())
{
@@ -233,10 +233,10 @@ double get_max_val(grid_type & 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)
+template <size_t Prop, typename prop_type, typename grid_type>
+prop_type get_min_val(grid_type & grid)
{
- double min_value = std::numeric_limits<double>::max();
+ prop_type min_value = std::numeric_limits<prop_type>::max();
auto dom = grid.getDomainIterator();
while(dom.isNext())
{
@@ -277,14 +277,14 @@ void init_sign_prop(grid_type & grid)
* @tparam gridtype Type of input grid.
* @param grid Grid, on which the magnitude of gradient should be computed.
*/
-template <size_t Vector_in, size_t Magnitude_out, typename gridtype>
+template <size_t Vector_in, size_t Magnitude_out, typename magnitude_type, typename gridtype>
void get_vector_magnitude(gridtype & grid)
{
grid.template ghost_get<Vector_in>();
auto dom = grid.getDomainGhostIterator();
while(dom.isNext())
{
- double sum = 0;
+ magnitude_type sum = 0;
auto key = dom.get();
for(size_t d = 0; d < gridtype::dims; d++)
{
diff --git a/src/level_set/redistancing_Sussman/RedistancingSussman.hpp b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
index 671c40325c4c04d5a9f015dd94244816c63d4f39..a95d0643dc7bbe9d65d14372b4eb8cf7f63a6430 100644
--- a/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
+++ b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
@@ -56,13 +56,14 @@
* between the iterations is considered and the
* redistancing finishes when the Conv_tol_change.value (or the max-iteration) is reached.
*/
+template <typename phi_type>
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-15; ///< Double variable that stores the tolerance value for the change at which the iterative
+ phi_type value = 1e-15; ///< Variable that stores the tolerance value for the change at which the iterative
///< redistancing process is considered as steady-state. (see also DistFromSol::change)
};
@@ -71,12 +72,13 @@ struct Conv_tol_change
* @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.
*/
+template <typename phi_type>
struct Conv_tol_residual
{
bool check = true; ///< 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: true.
- double value = 1e-3; ///< Double variable that stores the tolerance value for the residual at which the
+ phi_type value = 1e-3; ///< Variable that stores the tolerance value for the residual at which the
///< iterative redistancing process is considered as steady-state. (see also
///< DistFromSol::residual)
};
@@ -112,13 +114,14 @@ struct Conv_tol_residual
* @param save_temp_grid: If true, save the temporary grid as hdf5 that can be reloaded onto a grid
*/
+template <typename phi_type>
struct Redist_options
{
size_t min_iter = 1e5;
size_t max_iter = 1e12;
- Conv_tol_change convTolChange;
- Conv_tol_residual convTolResidual;
+ Conv_tol_change<phi_type> convTolChange;
+ Conv_tol_residual<phi_type> convTolResidual;
size_t interval_check_convergence = 100;
size_t width_NB_in_grid_points = 8;
@@ -135,13 +138,13 @@ struct Redist_options
*/
struct DistFromSol
{
- double change; ///< Double variable that contains the absolute value of the change of \a φ between the
+ auto change; ///< 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.
- double residual; ///< Double variable that contains the absolute value of how far the gradient magnitude of
+ auto residual; ///< 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.
@@ -154,7 +157,7 @@ struct DistFromSol
* @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>
+template <typename grid_in_type, typename phi_type>
class RedistancingSussman
{
public:
@@ -186,7 +189,7 @@ public:
* gradient of Phi_{n+1},
* sign of the original input Phi_0 (for the upwinding).
*/
- typedef aggregate<double, Point<grid_in_type::dims, double>, int>
+ typedef aggregate<phi_type, Point<grid_in_type::dims, phi_type>, int>
props_temp;
/** @brief Type definition for the temporary grid.
*/
@@ -235,7 +238,7 @@ public:
/** @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()
+ auto get_time_step()
{
/// This timestep is computed according to the grid spacing fulfilling the CFL condition.
return time_step;
@@ -246,12 +249,12 @@ public:
return final_iter;
}
- double get_finalChange()
+ auto get_finalChange()
{
return distFromSol.change;
}
- double get_finalResidual()
+ auto get_finalResidual()
{
return distFromSol.residual;
}
@@ -275,11 +278,11 @@ private:
int final_iter = 0; ///< Will be set to the final iteration when redistancing ends.
/// 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) * get_biggest_spacing(g_temp);
+ phi_type kappa = ceil(redistOptions.width_NB_in_grid_points / 2.0) * get_biggest_spacing(g_temp);
/**@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;
+ grid_type::stype time_step;
int order_upwind_gradient;
// Member functions
#ifdef SE_CLASS1
@@ -303,7 +306,7 @@ private:
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
+ phi_type 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_n_temp>(g_temp, min_value); // init. Phi_n_temp (incl. ghost) with min. Phi_0
copy_gridTogrid<Phi_0_in, Phi_n_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_n_temp
}
@@ -318,7 +321,7 @@ private:
* @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 sgn_phi_n)
+ phi_type get_phi_nplus1(phi_type phi_n, phi_type phi_n_magnOfGrad, grid_in_type::stype dt, phi_type sgn_phi_n)
{
return phi_n + dt * sgn_phi_n * (1 - phi_n_magnOfGrad);
}
@@ -335,9 +338,9 @@ private:
while (dom.isNext())
{
auto key = dom.get();
- const double phi_n = grid.template get<Phi_n_temp>(key);
- const double phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
- double epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
+ const phi_type phi_n = grid.template get<Phi_n_temp>(key);
+ const phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+ phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
grid.template get<Phi_n_temp>(key) = get_phi_nplus1(phi_n, phi_n_magnOfGrad, time_step,
smooth_S(phi_n, epsilon));
++dom;
@@ -350,7 +353,7 @@ private:
*
* @return True, if node lays within nb., false, if the distance to the interface is > kappa.
*/
- bool lays_inside_NB(double Phi)
+ bool lays_inside_NB(phi_type Phi)
{
return (abs(Phi) <= kappa);
}
@@ -362,8 +365,8 @@ private:
*/
void update_distFromSol(g_temp_type &grid)
{
- double max_residual = 0;
- double max_change = 0;
+ phi_type max_residual = 0;
+ phi_type max_change = 0;
int count = 0;
auto dom = grid.getDomainIterator();
while (dom.isNext())
@@ -372,9 +375,9 @@ private:
if (lays_inside_NB(grid.template get<Phi_n_temp>(key)))
{
count++;
- double phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
- double epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
- double phi_nplus1 = get_phi_nplus1(grid.template get<Phi_n_temp>(key), phi_n_magnOfGrad, time_step,
+ phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+ phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
+ phi_type phi_nplus1 = get_phi_nplus1(grid.template get<Phi_n_temp>(key), phi_n_magnOfGrad, time_step,
smooth_S(grid.template get<Phi_n_temp>(key), epsilon));
if (abs(phi_nplus1 - grid.template get<Phi_n_temp>(key)) > max_change)
diff --git a/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp b/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
index 7a75c9e4fa7f19afdb3a789e7c12576c54ea207f..c3d9516ebdd3942615fa5729cabf6ffb24bf8cb0 100644
--- a/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
+++ b/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
@@ -18,7 +18,9 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
BOOST_AUTO_TEST_CASE(RedistancingSussman_unit_sphere)
{
- const double EPSILON = std::numeric_limits<double>::epsilon();
+ 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;
@@ -31,19 +33,19 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
const size_t Error_grid = 3;
size_t N = 128;
- const double dt = 0.000165334; // CFL-condition for N=128
+ const space_type dt = 0.000165334; // CFL-condition for N=128
const size_t sz[grid_dim] = {N, N, N};
- const double radius = 1.0;
- const double box_lower = -2.0;
- const double box_upper = 2.0;
- Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
+ 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<double, double, double, double> props;
- typedef grid_dist_id<grid_dim, double, props > grid_in_type;
+ 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 double center[grid_dim] = {0.5*(box_upper+box_lower),
+ const space_type center[grid_dim] = {0.5*(box_upper+box_lower),
0.5*(box_upper+box_lower),
0.5*(box_upper+box_lower)};
init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
@@ -61,7 +63,8 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
redist_options.print_current_iterChangeResidual = true;
redist_options.print_steadyState_iter = true;
- RedistancingSussman<grid_in_type> redist_obj(g_dist, redist_options); // Instantiation of Sussman-redistancing class
+ RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options); // Instantiation of
+ // Sussman-redistancing class
std::cout << "New CFL timestep = " << redist_obj.get_time_step() << std::endl;
// redist_obj.set_user_time_step(dt);
// std::cout << "dt set to = " << dt << std::endl;
@@ -79,9 +82,9 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
// 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<double> props_nb;
- typedef vector_dist<grid_dim, double, props_nb> vd_type;
- Ghost<grid_dim, double> ghost_vd(0);
+ 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;