main.cpp

Go to the documentation of this file.

//
// Created by jstark on 2020-05-18.
//
 
 
// Include level-set-method related header files
#include "util/PathsAndFiles.hpp"
#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
#include "level_set/redistancing_Sussman/NarrowBand.hpp"
#include "Draw/DrawSphere.hpp"
 
int main(int argc, char* argv[])
{
    typedef double phi_type;
    //  initialize library
    openfpm_init(&argc, &argv);
    // Set current working directory, define output paths and create folders where output will be saved
    std::string cwd                     = get_cwd();
    const std::string path_output       = cwd + "/output_sphere";
    create_directory_if_not_exist(path_output);
    
    // grid dimension
    constexpr size_t grid_dim = 3;
    // Some indices for the grid / grid properties
    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 Phi_SDF_grid           = 1;
 
    // Here we create a 3D grid that stores two properties:
    // Prop1: store the initial Phi;
    // Prop2: here the re-initialized Phi (signed distance function) will be written to in the re-distancing step
    const size_t sz[grid_dim] = {128, 128, 128};
    Box<grid_dim, phi_type> box({0.0, 0.0, 0.0}, {10.0, 10.0, 10.0});
    Ghost<grid_dim, long int> ghost(0);
    typedef aggregate<phi_type, phi_type> props;
    typedef grid_dist_id<grid_dim, phi_type, props > grid_in_type;
    grid_in_type g_dist(sz, box, ghost);
    g_dist.setPropNames({"Phi_0", "Phi_SDF"});
    
    
    // Now we initialize the grid with a filled sphere. Outside the sphere, the value of Phi_0 will be -1, inside +1.
    phi_type radius = 1.0; // Radius of the sphere
    init_grid_with_sphere<Phi_0_grid>(g_dist, radius, 5, 5, 5); // Initialize sphere onto grid, centered at (5, 5, 5)
 
    g_dist.write(path_output + "/grid_initial_sphere_preRedistancing_radius" + std::to_string((int)radius) , FORMAT_BINARY); // Save the sphere as vtk file
 
 
 
    // Now we want to convert the initial Phi into a signed distance function (SDF) with magnitude of gradient = 1.
    // For the initial re-distancing we use the Sussman method. First of all, we can set some redistancing options.
    Redist_options<phi_type> redist_options;
    redist_options.min_iter                             = 1e3;
    redist_options.max_iter                             = 1e4;
    
    redist_options.convTolChange.value                  = 1e-7;
    redist_options.convTolChange.check                  = true;
    redist_options.convTolResidual.value                = 1e-6; // is ignored if convTolResidual.check = false;
    redist_options.convTolResidual.check                = false;
    
    redist_options.interval_check_convergence           = 1e3;
    redist_options.width_NB_in_grid_points              = 10;
    redist_options.print_current_iterChangeResidual     = true;
    redist_options.print_steadyState_iter               = true;
    redist_options.save_temp_grid                       = true;
    
    RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options);   // Instantiation of
    // Sussman-redistancing class
    //  std::cout << "dt = " << 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, Phi_SDF_grid>();
    
    g_dist.write(path_output + "/grid_sphere_postRedistancing", FORMAT_BINARY); // Save the result as vtk file
    g_dist.save(path_output + "/grid_sphere_postRedistancing" + ".bin"); // Save the result as hdf5 file that can be 
    // reloaded onto an openFPM grid
 
    //  Get narrow band: Place particles on interface (narrow band width e.g. 2 grid points on each side of the interface)
    size_t bc[grid_dim] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
    // Create an empty vector to which narrow-band particles will be added. You can choose, how many properties you want.
    // Minimum is 1 property, to which the Phi_SDF can be written
    // In this example we chose 3 properties. The 1st for the Phi_SDF, the 2nd for the gradient of phi and the 3rd for
    // the magnitude of the gradient
    typedef aggregate<phi_type, Point<grid_dim, phi_type>, phi_type> props_nb;
    typedef vector_dist<grid_dim, phi_type, props_nb> vd_type;
        Ghost<grid_dim, phi_type> ghost_vd(0);
        vd_type vd_narrow_band(0, box, bc, ghost_vd);
    vd_narrow_band.setPropNames({"Phi_SDF", "Phi_grad", "Phi_magnOfGrad"});
    
    size_t thickness_of_narrowBand_in_grid_points = 6;
    NarrowBand<grid_in_type, phi_type> narrowBand(g_dist, redist_options.width_NB_in_grid_points); // Instantiation of
    // NarrowBand class
 
    
    // Some indices
    const size_t Phi_SDF_vd             = 0;
    const size_t Phi_grad_vd            = 1;
    const size_t Phi_magnOfGrad_vd      = 2;
    
    // Get the narrow band. You can decide, if you only want the Phi_SDF saved to your particles or
    // if you also want the gradients or gradients and magnitude of gradient.
    // The function will know depending on how many property-indices you provide in the <> brackets.
    // First property-index must always be the index of the SDF on the grid!
    // E.g.: The following line would only write only the Phi_SDF from the grid to your narrow-band particles
    // narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd>(g_dist, vd_narrow_band);
    // Whereas this will give you the gradients and magnOfGrad of Phi as well:
    narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd, Phi_grad_vd, Phi_magnOfGrad_vd>(g_dist, vd_narrow_band);
    
    vd_narrow_band.write(path_output + "/vd_narrow_band_sphere", FORMAT_BINARY); // Save particles as vtk file
    vd_narrow_band.save(path_output + "/vd_narrow_band_sphere.bin"); // Save particles as hdf5 file -> can be reloaded as particles
    
    openfpm_finalize(); // Finalize openFPM library
    return 0;
}