Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Sbalzarini Lab
S
Software
P
Parallel Computing
OpenFPM
openfpm_numerics
Commits
c16a2476
Commit
c16a2476
authored
Dec 30, 2021
by
jstark
Browse files
Generalizing types in Sussman redistancing and corresponding HelpFunctionsForGrid.
parent
d0d7fe8c
Changes
4
Hide whitespace changes
Inline
Side-by-side
src/CMakeLists.txt
View file @
c16a2476
...
...
@@ -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
)
...
...
src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
View file @
c16a2476
...
...
@@ -26,9 +26,9 @@
* @return Time step.
*/
template
<
typename
grid_type
>
doubl
e
get_time_step_CFL
(
grid_type
&
grid
)
typename
grid_type
::
styp
e
get_time_step_CFL
(
grid_type
&
grid
)
{
doubl
e
sum
=
0
;
typename
grid_type
::
styp
e
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>
doubl
e get_time_step_CFL(grid_type & grid,
double u [grid_type::dims], doubl
e C)
typename grid_type::styp
e get_time_step_CFL(grid_type & grid,
typename grid_type::stype u [grid_type::dims], typename grid_type::styp
e C)
{
doubl
e sum = 0;
typename grid_type::styp
e 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>
doubl
e get_time_step_CFL(grid_type & grid,
double u, doubl
e C)
typename grid_type::styp
e get_time_step_CFL(grid_type & grid,
typename grid_type::stype u, typename grid_type::styp
e C)
{
doubl
e sum = 0;
typename grid_type::styp
e 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
>
doubl
e
get_biggest_spacing
(
grid_type
&
grid
)
typename
grid_type
::
styp
e
get_biggest_spacing
(
grid_type
&
grid
)
{
doubl
e
h_max
=
0
;
typename
grid_type
::
styp
e
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
>
doubl
e
get_smallest_spacing
(
grid_type
&
grid
)
typename
grid_type
::
styp
e
get_smallest_spacing
(
grid_type
&
grid
)
{
doubl
e
spacing
[
grid_type
::
dims
];
typename
grid_type
::
styp
e
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
>
doubl
e
average_difference
(
grid_type
&
grid
)
template
<
size_t
Prop1
,
size_t
Prop2
,
typename
prop_type
,
typename
grid_type
>
prop_typ
e
average_difference
(
grid_type
&
grid
)
{
doubl
e
total_diff
=
0
;
prop_typ
e
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
>
doubl
e
get_max_val
(
grid_type
&
grid
)
template
<
size_t
Prop
,
typename
prop_type
,
typename
grid_type
>
prop_typ
e
get_max_val
(
grid_type
&
grid
)
{
doubl
e
max_value
=
std
::
numeric_limits
<
doubl
e
>::
lowest
();
prop_typ
e
max_value
=
std
::
numeric_limits
<
prop_typ
e
>::
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
>
doubl
e
get_min_val
(
grid_type
&
grid
)
template
<
size_t
Prop
,
typename
prop_type
,
typename
grid_type
>
prop_typ
e
get_min_val
(
grid_type
&
grid
)
{
doubl
e
min_value
=
std
::
numeric_limits
<
doubl
e
>::
max
();
prop_typ
e
min_value
=
std
::
numeric_limits
<
prop_typ
e
>::
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
())
{
doubl
e
sum
=
0
;
magnitude_typ
e
sum
=
0
;
auto
key
=
dom
.
get
();
for
(
size_t
d
=
0
;
d
<
gridtype
::
dims
;
d
++
)
{
...
...
src/level_set/redistancing_Sussman/RedistancingSussman.hpp
View file @
c16a2476
...
...
@@ -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.
doubl
e
value
=
1e-15
;
///<
Double v
ariable that stores the tolerance value for the change at which the iterative
phi_typ
e
value
=
1e-15
;
///<
V
ariable 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.
doubl
e
value
=
1e-3
;
///<
Double v
ariable that stores the tolerance value for the residual at which the
phi_typ
e
value
=
1e-3
;
///<
V
ariable 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 v
ariable that contains the absolute value of the change of \a φ between the
auto
change
;
///<
V
ariable 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 v
ariable that contains the absolute value of how far the gradient magnitude of
auto
residual
;
///<
V
ariable 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
<
doubl
e
,
Point
<
grid_in_type
::
dims
,
doubl
e
>
,
int
>
typedef
aggregate
<
phi_typ
e
,
Point
<
grid_in_type
::
dims
,
phi_typ
e
>
,
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.
doubl
e
kappa
=
ceil
(
redistOptions
.
width_NB_in_grid_points
/
2.0
)
*
get_biggest_spacing
(
g_temp
);
phi_typ
e
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()
*/
doubl
e
time_step
;
grid_type
::
styp
e
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
()
{
doubl
e
min_value
=
get_min_val
<
Phi_0_in
>
(
r_grid_in
);
// get minimum Phi_0 value on the input grid
phi_typ
e
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.
*
*/
doubl
e
get_phi_nplus1
(
doubl
e
phi_n
,
doubl
e
phi_n_magnOfGrad
,
double
dt
,
doubl
e
sgn_phi_n
)
phi_typ
e
get_phi_nplus1
(
phi_typ
e
phi_n
,
phi_typ
e
phi_n_magnOfGrad
,
grid_in_type
::
stype
dt
,
phi_typ
e
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
doubl
e
phi_n
=
grid
.
template
get
<
Phi_n_temp
>(
key
);
const
doubl
e
phi_n_magnOfGrad
=
grid
.
template
get
<
Phi_grad_temp
>(
key
).
norm
();
doubl
e
epsilon
=
phi_n_magnOfGrad
*
grid
.
getSpacing
()[
0
];
const
phi_typ
e
phi_n
=
grid
.
template
get
<
Phi_n_temp
>(
key
);
const
phi_typ
e
phi_n_magnOfGrad
=
grid
.
template
get
<
Phi_grad_temp
>(
key
).
norm
();
phi_typ
e
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
(
doubl
e
Phi
)
bool
lays_inside_NB
(
phi_typ
e
Phi
)
{
return
(
abs
(
Phi
)
<=
kappa
);
}
...
...
@@ -362,8 +365,8 @@ private:
*/
void
update_distFromSol
(
g_temp_type
&
grid
)
{
doubl
e
max_residual
=
0
;
doubl
e
max_change
=
0
;
phi_typ
e
max_residual
=
0
;
phi_typ
e
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
++
;
doubl
e
phi_n_magnOfGrad
=
grid
.
template
get
<
Phi_grad_temp
>(
key
).
norm
();
doubl
e
epsilon
=
phi_n_magnOfGrad
*
grid
.
getSpacing
()[
0
];
doubl
e
phi_nplus1
=
get_phi_nplus1
(
grid
.
template
get
<
Phi_n_temp
>(
key
),
phi_n_magnOfGrad
,
time_step
,
phi_typ
e
phi_n_magnOfGrad
=
grid
.
template
get
<
Phi_grad_temp
>(
key
).
norm
();
phi_typ
e
epsilon
=
phi_n_magnOfGrad
*
grid
.
getSpacing
()[
0
];
phi_typ
e
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
)
...
...
src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
View file @
c16a2476
...
...
@@ -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
doubl
e
dt
=
0.000165334
;
// CFL-condition for N=128
const
space_typ
e
dt
=
0.000165334
;
// CFL-condition for N=128
const
size_t
sz
[
grid_dim
]
=
{
N
,
N
,
N
};
const
doubl
e
radius
=
1.0
;
const
doubl
e
box_lower
=
-
2.0
;
const
doubl
e
box_upper
=
2.0
;
Box
<
grid_dim
,
doubl
e
>
box
({
box_lower
,
box_lower
,
box_lower
},
{
box_upper
,
box_upper
,
box_upper
});
const
space_typ
e
radius
=
1.0
;
const
space_typ
e
box_lower
=
-
2.0
;
const
space_typ
e
box_upper
=
2.0
;
Box
<
grid_dim
,
space_typ
e
>
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
,
doubl
e
>
props
;
typedef
grid_dist_id
<
grid_dim
,
doubl
e
,
props
>
grid_in_type
;
typedef
aggregate
<
phi_type
,
phi_type
,
phi_type
,
phi_typ
e
>
props
;
typedef
grid_dist_id
<
grid_dim
,
space_typ
e
,
props
>
grid_in_type
;
grid_in_type
g_dist
(
sz
,
box
,
ghost
);
g_dist
.
setPropNames
({
"Phi_0"
,
"SDF_sussman"
,
"SDF_exact"
,
"Relative error"
});
const
doubl
e
center
[
grid_dim
]
=
{
0.5
*
(
box_upper
+
box_lower
),
const
space_typ
e
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
<
doubl
e
>
props_nb
;
typedef
vector_dist
<
grid_dim
,
doubl
e
,
props_nb
>
vd_type
;
Ghost
<
grid_dim
,
doubl
e
>
ghost_vd
(
0
);
typedef
aggregate
<
phi_typ
e
>
props_nb
;
typedef
vector_dist
<
grid_dim
,
space_typ
e
,
props_nb
>
vd_type
;
Ghost
<
grid_dim
,
space_typ
e
>
ghost_vd
(
0
);
vd_type
vd_narrow_band
(
0
,
box
,
bc
,
ghost_vd
);
vd_narrow_band
.
setPropNames
({
"error"
});
size_t
narrow_band_width
=
8
;
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment