NarrowBand.hpp

Go to the documentation of this file.

//
// Created by jstark on 2020-05-16.
//
#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 "HelpFunctions.hpp"
#include "HelpFunctionsForGrid.hpp"
#include "FiniteDifference/Upwind_gradient.hpp"

template <typename grid_in_type, typename phi_type=double>
class NarrowBand
{
public:
    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);
    }
    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);
    }
    NarrowBand(const grid_in_type & grid_in,
               float 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);
    }
    template<typename width_type>
    NarrowBand(const grid_in_type & grid_in,
               width_type width_outside,     // physical width of nb inside of object -> Phi > 0
               width_type 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);
    }
    typedef aggregate<phi_type, phi_type[grid_in_type::dims], int> 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;
    
    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)
    {
        auto dom = grid.getDomainIterator();
        while(dom.isNext())
        {
            auto key = dom.get();
            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] = grid.getPos(key)[d];
                }
                vd.template getLastProp<Phi_SDF_vd>() = grid.template get<Phi_SDF_grid>(key);
            }
            ++dom;
        }
        vd.map();
    }
    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)
    {
        initialize_temporary_grid<Phi_SDF_grid>(grid);
        auto dom = grid.getDomainIterator();
        while(dom.isNext())
        {
            auto key = dom.get();
            if(within_narrow_band(grid.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] = grid.getPos(key)[d];
                    vd.template getLastProp<Phi_grad>()[d] = g_temp.template get<Phi_grad_temp>(key)[d];
                }
                vd.template getLastProp<Phi_SDF_vd>() = grid.template get<Phi_SDF_temp>(key);
            }
            ++dom;
        }
        vd.map();
    }
    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)
    {
        initialize_temporary_grid<Phi_SDF_grid>(grid);
        auto dom = grid.getDomainIterator();
        while(dom.isNext())
        {
            auto key = dom.get();
            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] = grid.getPos(key)[d];
                    vd.template getLastProp<Phi_grad>()[d] = g_temp.template get<Phi_grad_temp>(key)[d];
                }
                vd.template getLastProp<Phi_SDF_vd>()      = g_temp.template get<Phi_SDF_temp>(key);
                vd.template getLastProp<Phi_magnOfGrad>()  = get_vector_magnitude<Phi_grad_temp>(g_temp, key);
            }
            ++dom;
        }
        vd.map();   }
    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)
    {
        auto dom = grid.getDomainIterator();
        while(dom.isNext())
        {
            auto key = dom.get();
            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] = grid.getPos(key)[d];
                }
                vd.template getLastProp<Prop1_vd>() = grid.template get<Prop1_grid>(key);
            }
            ++dom;
        }
        vd.map();
    }
    template <size_t Phi_SDF_grid, size_t Index1Grid, size_t Index2Grid, size_t Index3Grid, size_t Index1Vd, size_t
    Index2Vd, size_t Index3Vd, typename grid_type, typename vector_type>
    void get_narrow_band_copy_three_scalar_properties(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<Index1Vd>() = grid.template get<Index1Grid>(key);
                    vd.template getLastProp<Index2Vd>() = grid.template get<Index2Grid>(key);
                    vd.template getLastProp<Index3Vd>() = grid.template get<Index3Grid>(key);
                }
                ++dom;
            }
            vd.map();
    }
private:
    //  Some indices for better readability
    static const size_t Phi_SDF_temp        = 0; 
    static const size_t Phi_grad_temp       = 1; 
    static const size_t Phi_sign_temp       = 2; 
    
    phi_type b_low;   
    phi_type b_up;    
    
    void set_bounds(const size_t thickness, const grid_in_type & grid_in)
    {
        const phi_type EPSILON = std::numeric_limits<phi_type>::epsilon();
        b_low   = - ceil((phi_type)thickness / 2.0) * get_biggest_spacing(grid_in) - EPSILON;
        b_up    =   ceil((phi_type)thickness / 2.0) * get_biggest_spacing(grid_in) + EPSILON;
    }
    template<typename thickness_type>
    void set_bounds(const thickness_type thickness, const grid_in_type & grid_in)
    {
        const thickness_type EPSILON = std::numeric_limits<thickness_type>::epsilon();
        b_low   = - thickness / 2.0 - EPSILON;
        b_up    =   thickness / 2.0 + EPSILON;
    }
    template<typename width_type>
    void set_bounds(width_type lower_bound, width_type upper_bound, const grid_in_type & grid_in)
    {
        const width_type EPSILON = std::numeric_limits<width_type>::epsilon();
        b_low   =   lower_bound - EPSILON;
        b_up    =   upper_bound + EPSILON;
    }
    
    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, 1, true);   // Get initial gradients
    }
    template<typename T>
    bool within_narrow_band(T Phi)
    {
        return (Phi >= b_low && Phi <= b_up);
    }
};


#endif //REDISTANCING_SUSSMAN_NARROWBAND_HPP