#ifndef RCFD_USER
#define RCFD_USER

#include "udf.h"
#include "rCFD_library.h"

/* (C)  2023
    Stefan Pirker
    Particulate Flow Modelling
    Johannes Kepler University, Linz, Austria
    www.particulate-flow.at
*/


    /* user Names */
    /*************************************************************************************/

#if 1

#if RP_NODE

    #define     _diam_23    (2.0 * 0.001142)
    
    #define     _k_diam     0.2
    
    #define     _number_of_particle_fractions   6
    
    #define     _diam_0     (2. * _diam_23)
    #define     _diam_1     (3./2. * _diam_23)
    #define     _diam_2     (5./4. * _diam_23)
    #define     _diam_3     _diam_23
    #define     _diam_4     (7./8. * _diam_23)
    #define     _diam_5     (3./4. * _diam_23)
    
    #define     _mf_diam_0  0.0
    #define     _mf_diam_1  0.1
    #define     _mf_diam_2  0.15
    #define     _mf_diam_3  0.445
    #define     _mf_diam_4  0.205
    #define     _mf_diam_5  0.1

    #define     loop_diam           for(i_diam = 0; i_diam < _number_of_particle_fractions; i_diam ++)
    #define     loop_int_cells      for(i_cell = 0; i_cell < int_cell_count; i_cell++)
    
    static short    rCFD_monitors_initiated = 0;
    
#endif

    /* names of phases */
    enum{
        ph_gas,
        ph_solid,
        number_of_phase_names
    };

    /* names of data - phase_0 */
    enum{
        number_of_phase0_data_names
    };

    /* names of data - phase_1 */
    enum{
        c_r0,
        c_r1,
        c_r2,
        c_r3,
        c_r4,
        c_r5,
        
        number_of_phase1_data_names
    };

    enum{
        c_vf_diam_0,
        c_vf_diam_1,
        c_vf_diam_2,
        c_vf_diam_3,
        c_vf_diam_4,
        c_vf_diam_5,
        
        c_vof,

        c_vf_diam_0_ref,
        c_vf_diam_1_ref,
        c_vf_diam_2_ref,
        c_vf_diam_3_ref,
        c_vf_diam_4_ref,
        c_vf_diam_5_ref,

        c_vof_ref,
        
        c_vf_diam_0_av,
        c_vf_diam_1_av,
        c_vf_diam_2_av,
        c_vf_diam_3_av,
        c_vf_diam_4_av,
        c_vf_diam_5_av,

        c_vf_diam_0_ref_av,
        c_vf_diam_1_ref_av,
        c_vf_diam_2_ref_av,
        c_vf_diam_3_ref_av,
        c_vf_diam_4_ref_av,
        c_vf_diam_5_ref_av,

        number_of_UDMI_names
    };

    #define     _number_of_bins     25
    
    #define     _y_min      0.00
    #define     _y_max      0.25
    
    #define     _y_bottom   0.06

    
#if RP_NODE
    
    static  double  part_diam[6];       // set by rCFD_user_init_Data

    static  double  part_mf_diam[6];    // set by rCFD_user_init_Data
    
    int _growth_zone(double *x)
    {
        if((x[0] > 0.09) && (x[0] < 0.11) && (x[1] > 0.06)){
            
            return 1;
        }
        else{
            
            return 0;
        }
    }

    int _inflow_zone(double *x)
    {
        if((x[0] < 0.025) && (x[1] < 0.06)){
            
            return 1;
        }
        else{
            
            return 0;
        }
    }

    int _outflow_zone(double *x)
    {
        if((x[0] > 0.175) && (x[1] < 0.06)){
            
            return 1;
        }
        else{
            
            return 0;
        }
    }
    
    #define     _global_growth_rate    0.5      // 5.0e-1    // (kg/s)
    
    #define     _cross_fract_coef      0.1
    
    #define     _global_throughflow    0.25     // 2.5e-1    // (kg/s)

#endif
    
#endif

    /* user set Dicts */
    /*************************************************************************************/

#if 1

    void user_set_solver_dict(void *handle)
    {
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "report_priority", 1);
        
        rcfd_set_dict_property_r(handle, "Solver_Dict", -1, -1, -1, "global_time_step", (2.0 * 5.0e-3));

        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "number_of_frames", (300 / 2));

        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "number_of_phases", 2);

        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "face_swap_diffusion_on", 0);
        
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "data_drifting_on", 1);
        
        rcfd_set_dict_property_r(handle, "Solver_Dict", -1, -1, -1, "vof_cutoff_limit", 1.0e-6);
    }

    void user_set_phase_dict(void *handle)
    {
#if RP_NODE

        rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_gas, -1, "number_of_data", 0);

        rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "number_of_data", 6);
        
        rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "density", 1450.0);
        
        rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "vof_max", 0.65);

        rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "d23_diam", _diam_23);
        
        rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "hindered_drift_on", 0);   

        rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "local_control_conc_sum_on", 1);   

        rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "local_drift_on", 0);  
        
#endif
    }

    void user_set_rec_dict(void *handle)
    {
        rcfd_set_dict_property_s(handle, "Rec_Dict", -1, -1, -1, "format", "quarter_jumps_format");
    }

    void user_set_data_dict(void *handle)
    {
#if RP_NODE
        int i_diam;
        
        loop_diam{
            
            rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, i_diam, "type", "concentration_data");
            
            rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, i_diam, "drifting_data", 1);
        }

        double  drift_velocity[3];
        
        drift_velocity[0] = 0.0;
        
        drift_velocity[2] = 0.0;
        
        // c_r0 = LLarge particles
        {
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r0, "d1_diam", _diam_0);
            
            drift_velocity[1] = - _k_diam * (_diam_0 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r0, "drift_velocity", drift_velocity);
        }

        // c_r1 = Large particles
        {
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r1, "d1_diam", _diam_1);
            
            drift_velocity[1] = - _k_diam * (_diam_1 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r1, "drift_velocity", drift_velocity);
        }
        
        // c_r2 = large particles
        {           
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r2, "d1_diam", _diam_2);
            
            drift_velocity[1] = - _k_diam * (_diam_2 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r2, "drift_velocity", drift_velocity);
        }

        // c_r3 = d23 particles
        {
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r3, "d1_diam", _diam_3);
            
            drift_velocity[1] = - _k_diam * (_diam_3 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r3, "drift_velocity", drift_velocity);
        }

        // c_r4 = small particles
        {
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r4, "d1_diam", _diam_4);
            
            drift_velocity[1] = - _k_diam * (_diam_4 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r4, "drift_velocity", drift_velocity);
        }

        // c_r5 = Small particles
        {
            rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r5, "d1_diam", _diam_5);
            
            drift_velocity[1] = - _k_diam * (_diam_5 - _diam_23);
            
            rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r5, "drift_velocity", drift_velocity);
        }
        
#endif
    }

    void user_set_cell_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        fluent_set_cell_dict(handle, i_layer);
#endif
    }

    void user_set_face_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        fluent_set_face_dict(handle, i_layer);
#endif
    }   

#endif

    /* user init */
    /*************************************************************************************/

#if 1

    void user_init_data(void *handle, const short i_layer)
    {
#if RP_NODE

        int     i_cell;

        int     int_cell_count  = rcfd_get_int_cell_count(handle, i_layer);
        
        double  **solid_data = rcfd_get_cell_data(handle, i_layer, ph_solid);       

        loop_int_cells{

            solid_data[i_cell][c_r0] = _mf_diam_0;            
            solid_data[i_cell][c_r1] = _mf_diam_1;            
            solid_data[i_cell][c_r2] = _mf_diam_2;
            solid_data[i_cell][c_r3] = _mf_diam_3;            
            solid_data[i_cell][c_r4] = _mf_diam_4;
            solid_data[i_cell][c_r5] = _mf_diam_5;            
        }
        
        part_diam[0] = _diam_0;
        part_diam[1] = _diam_1;
        part_diam[2] = _diam_2;
        part_diam[3] = _diam_3;
        part_diam[4] = _diam_4;
        part_diam[5] = _diam_5;
        
        part_mf_diam[0] = _mf_diam_0;
        part_mf_diam[1] = _mf_diam_1;
        part_mf_diam[2] = _mf_diam_2;
        part_mf_diam[3] = _mf_diam_3;
        part_mf_diam[4] = _mf_diam_4;
        part_mf_diam[5] = _mf_diam_5;

        Message0("\n\n...user_init_Data\n");
#endif
    }

#endif

    /* user accesses Data */
    /*************************************************************************************/

#if 1

    void user_access_data_before_phase_loop(void *handle)
    {
#if RP_NODE
  
#if 1   // local vars
        
        int     i_layer, i_cell, i_diam, i_frame;
        
        int     int_cell_count; 
        
        double  local_part_vol, local_fract_part_vol, local_fract_part_area;

        double  global_time_step, vof_cutoff_limit, solid_density;
        
        double  per_node_part_area, global_growth_area, cross_fract_exchange;
        
        double  local_mass, local_fract_inflow, sum_of_conc, per_node_inlet_part_mass, global_inlet_part_mass;
        
        double  per_node_outlet_part_mass, global_outlet_part_mass;
        
        double  *per_node_outlet_part_fract_mass = NULL, *global_outlet_part_fract_mass = NULL, *average_part_fraction = NULL;
        
        double  *balance_source = NULL, local_balance_source;
        
        double  *cell_volume = NULL;

        double  **x = NULL, **solid_data = NULL, **vof = NULL;
        
#endif      
        
        // set local vars, access global vars
        {
            i_layer = 0;
            
            i_frame = rcfd_get_rec_global_frame(handle);

            solid_data = rcfd_get_cell_data(handle, i_layer, ph_solid);
            
            solid_density = rcfd_get_phase_density(handle, ph_solid);
    
            int_cell_count  = rcfd_get_int_cell_count(handle, i_layer);

            x = rcfd_get_cell_center(handle, i_layer); // [i_cell][i_dim]
                        
            cell_volume = rcfd_get_cell_volume(handle, i_layer);

            vof = rcfd_get_cell_vof(handle, i_layer, i_frame);  // [i_cell][i_phase]
            
            vof_cutoff_limit = rcfd_get_vof_cutoff_limit(handle);
            
            global_time_step = rcfd_get_timestep(handle, i_layer);
        }
        
        // calc. global_growth_area of particles in growth zone
        {
            per_node_part_area = 0.0;
            
            loop_int_cells{
                
                if(_growth_zone(x[i_cell])){
                    
                    local_part_vol = vof[i_cell][ph_solid] * cell_volume[i_cell];
                    
                    loop_diam{
                        
                        local_fract_part_vol = solid_data[i_cell][i_diam] * local_part_vol;
                        
                        local_fract_part_area = 3.0/2.0 * local_fract_part_vol / part_diam[i_diam];
                        
                        per_node_part_area += local_fract_part_area;
                    }
                }
            }
            
            global_growth_area = PRF_GRSUM1(per_node_part_area);
        }
        
        // adapt fractional concentrations based on total growth rate
        {
            
            balance_source = (double*)malloc(_number_of_particle_fractions * sizeof(double));
            
            loop_diam{
                
                balance_source[i_diam] = 0.0;
            }
            
            loop_int_cells{
                
                if(_growth_zone(x[i_cell])){
                    
                    if(vof[i_cell][ph_solid] > vof_cutoff_limit){
                        
                        local_part_vol = vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                        
                            local_fract_part_vol = solid_data[i_cell][i_diam] * local_part_vol;
                            
                            local_fract_part_area = 3.0/2.0 * local_fract_part_vol / part_diam[i_diam];
                            
                            if( i_diam > 0){    // i_diam_0 is largest fraction ...
                                
                                cross_fract_exchange = _cross_fract_coef * local_fract_part_area / global_growth_area *
                                
                                    _global_growth_rate * global_time_step / solid_density;  // (m3)
                                    
                                if(cross_fract_exchange < solid_data[i_cell][i_diam] * local_part_vol){
                                    
                                    solid_data[i_cell][i_diam] = (solid_data[i_cell][i_diam] * local_part_vol - cross_fract_exchange) / local_part_vol;
                                }
                                else{
                                    
                                    cross_fract_exchange = solid_data[i_cell][i_diam] * local_part_vol;
                                    
                                    solid_data[i_cell][i_diam] = 0.0;
                                }
                                
                                solid_data[i_cell][(i_diam - 1)] = (solid_data[i_cell][(i_diam - 1)] * local_part_vol + cross_fract_exchange) / local_part_vol;
                                
                                local_balance_source = cross_fract_exchange * solid_density;    // (kg)
                                
                                balance_source[i_diam] -=  local_balance_source;                //  smaller particles loose mass ...
                                
                                balance_source[(i_diam - 1)] += local_balance_source;   
                            }
                        }
                    }
                }
            }

            rcfd_modify_balance_mass_source_each(handle, ph_solid, balance_source);         
        }
        
        // calc global_part_mass[i_diam] (i.e. mass of particles of fraction i_diam within inlet zone)
        {
            
            per_node_inlet_part_mass = 0.0;
            
            loop_int_cells{
                
                if(_inflow_zone(x[i_cell])){
                
                    per_node_inlet_part_mass += vof[i_cell][ph_solid] * cell_volume[i_cell] * solid_density;
                }
            }
            
            global_inlet_part_mass = PRF_GRSUM1(per_node_inlet_part_mass);
        }
        
        // adjust solid_data and balance according to inflow
        {
            loop_diam{
                
                balance_source[i_diam] = 0.0;
            }
            
            loop_int_cells{
                
                if(_inflow_zone(x[i_cell])){         
                
                    if(vof[i_cell][ph_solid] > vof_cutoff_limit){
                        
                        local_mass = vof[i_cell][ph_solid] * cell_volume[i_cell] * solid_density;
                        
                        loop_diam{

                            if(global_inlet_part_mass > 0.0){
                                
                                local_fract_inflow = local_mass / global_inlet_part_mass *
                                
                                    _global_throughflow * part_mf_diam[i_diam] * global_time_step;
                            
                                solid_data[i_cell][i_diam] = (solid_data[i_cell][i_diam] * local_mass + local_fract_inflow) / local_mass;
                            
                                balance_source[i_diam] += local_fract_inflow;
                            }
                        }
                    }
                }
            }
            
            rcfd_modify_balance_mass_source_each(handle, ph_solid, balance_source);                     
        }

        // adjust solid_data such that their sum is one
        {
            loop_int_cells{
                
                if(_inflow_zone(x[i_cell])){ 
                
                    sum_of_conc = 0.0;
                    
                    loop_diam{
                            
                        sum_of_conc += solid_data[i_cell][i_diam];
                    }
                    
                    if( sum_of_conc <= 0.0) sum_of_conc = 1.0;  // don't change anything
                    
                    loop_diam{
    
                        solid_data[i_cell][i_diam] /= sum_of_conc;
                    }
                }
            }
        }
        
        // account for outflow of particles 
        {
            // initializations
            {
                per_node_outlet_part_fract_mass =   (double*)malloc(_number_of_particle_fractions * sizeof(double));
                
                global_outlet_part_fract_mass =     (double*)malloc(_number_of_particle_fractions * sizeof(double));
                
                average_part_fraction =             (double*)malloc(_number_of_particle_fractions * sizeof(double));
                
                loop_diam{
                    
                    per_node_outlet_part_fract_mass[i_diam] =   0.0;
                    
                    global_outlet_part_fract_mass[i_diam] =     0.0;
                    
                    average_part_fraction[i_diam] =             0.0;
                }
                
                per_node_outlet_part_mass = 0.0;
            }
            
            // calc global_outlet_part_fract_mass[i_diam]
            {
                loop_int_cells{
                    
                    if(_outflow_zone(x[i_cell])){
                    
                        per_node_outlet_part_mass += vof[i_cell][ph_solid] * cell_volume[i_cell] * solid_density;
                        
                        loop_diam{
                            
                            per_node_outlet_part_fract_mass[i_diam] += solid_data[i_cell][i_diam] * vof[i_cell][ph_solid] * cell_volume[i_cell] * solid_density;
                        }
                    }
                }
                
                global_outlet_part_mass = PRF_GRSUM1(per_node_outlet_part_mass);
                
                loop_diam{
                    
                    global_outlet_part_fract_mass[i_diam] = PRF_GRSUM1(per_node_outlet_part_fract_mass[i_diam]);
                }
            }
            
            // adapt Balances
            if(global_outlet_part_mass > 0.0){
                
                loop_diam{
                    
                    average_part_fraction[i_diam] = global_outlet_part_fract_mass[i_diam] / global_outlet_part_mass;
                
                    balance_source[i_diam] = -average_part_fraction[i_diam] * per_node_outlet_part_mass / global_outlet_part_mass * 
                    
                        _global_throughflow * global_time_step;    // (kg)
                }
                
                rcfd_modify_balance_mass_source_each(handle, ph_solid, balance_source); 
            }
        
        }

        // free local vars
        {
            
            free(balance_source);
            
            free(per_node_outlet_part_fract_mass);

            free(global_outlet_part_fract_mass);

            free(average_part_fraction);            
        }
        
#endif
    }

    void user_post(void *handle)
    {
#if RP_NODE

#if 1   /* local vars */
        int     i_layer, i_cell, i_frame, i_diam;
        
        int     int_cell_count;

        double  vol_bottom, vol_top, global_time;
        
        double  *vol_bottom_diam = NULL, *vol_top_diam = NULL;
        
        double  *mf_bottom_diam = NULL, *mf_top_diam = NULL;
        
        double  *cell_volume = NULL;
        
        double  **x = NULL, **solid_data = NULL, **vof = NULL;

        Domain  *d = Get_Domain(1);

        Thread  *t = NULL, *t_mix = NULL;
#endif

        // access vars, set local vars
        {
            i_layer = 0;
            
            int_cell_count  = rcfd_get_int_cell_count(handle, i_layer);
            
            i_frame = rcfd_get_rec_global_frame(handle);

            solid_data = rcfd_get_cell_data(handle, i_layer, ph_solid);
            
            x = rcfd_get_cell_center(handle, i_layer); // [i_cell][i_dim]
            
            vof = rcfd_get_cell_vof(handle, i_layer, i_frame);  // [i_cell][i_phase]
            
            cell_volume = rcfd_get_cell_volume(handle, i_layer);

            global_time = rcfd_get_time(handle);
            
            // get Fluent t_mix
            thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){

                t_mix = t;
            }}
        }   
        
        // set UDMI(vf_diam)
        {
            loop_int_cells{
                
                loop_diam{
                    
                    C_UDMI(i_cell, t_mix, (c_vf_diam_0 + i_diam)) = solid_data[i_cell][i_diam] * vof[i_cell][ph_solid];
                }
                
                C_UDMI(i_cell, t_mix, (c_vf_diam_0 + 6)) = vof[i_cell][ph_solid];
            }
        }
        
        // write bottom/top psd
        {
            // bottom / top masses (rCFD, ref)
            {
                // initializations
                {
                    vol_bottom      = 0.0;              
                    vol_top         = 0.0;
                    
                    vol_bottom_diam =       (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    vol_top_diam =          (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] =       0.0;
                        vol_top_diam[i_diam] =          0.0;
                    }
                }
                
                loop_int_cells{
                
                    if(x[i_cell][1] < _y_bottom){
                        
                        // rCFD
                        vol_bottom += vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_bottom_diam[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_0 + i_diam)) * cell_volume[i_cell];
                        }
                    }
                    else{
                        
                        // rCFD
                        vol_top += vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_top_diam[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_0 + i_diam)) * cell_volume[i_cell];
                        }
                    }
                }
            
                // global reduction
                {
                    vol_bottom =        PRF_GRSUM1(vol_bottom);
                    vol_top =           PRF_GRSUM1(vol_top);
                     
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] =       PRF_GRSUM1(vol_bottom_diam[i_diam]);
                        vol_top_diam[i_diam] =          PRF_GRSUM1(vol_top_diam[i_diam]);
                    }   
                }
                
                // mf_diam top/bottom
                {
                    // init mf's
                    {
                        mf_bottom_diam =        (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        mf_top_diam =           (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        
                        loop_diam{
                            
                            mf_bottom_diam[i_diam] =        0.0;
                            mf_top_diam[i_diam] =           0.0;
                        }
                    }
                    
                    if(vol_bottom > 0.0){
                        
                        loop_diam{
                            
                            mf_bottom_diam[i_diam] = vol_bottom_diam[i_diam] / vol_bottom;
                        }
                    }
                    
                    if(vol_top > 0.0){
                        
                        loop_diam{
                            
                            mf_top_diam[i_diam] = vol_top_diam[i_diam] / vol_top;
                        }
                    }
                }
            }
            
            // write monitor file
            if(myid == 0){

                FILE    *f_out = NULL;

                if(rCFD_monitors_initiated == 0){

                    f_out = fopen("./rCFD_monitor_poly.out", "w");

                    rCFD_monitors_initiated = 1;
                }
                else{
                    f_out = fopen("./rCFD_monitor_poly.out", "a");
                }

                if(f_out == NULL){

                    Message0("\nERROR: Could not open monitor file\n");

                    return;
                }

                fprintf(f_out, "%e ", global_time);
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_bottom_diam[i_diam]);
                }
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_top_diam[i_diam]);
                }
                
                fprintf(f_out, "\n");

                fclose(f_out);
            }
        
        }
        
        // free local vars
        {
    
            free(vol_bottom_diam);
            free(vol_top_diam);
    
            free(mf_bottom_diam);
            free(mf_top_diam);
         }

        Message0("\n\n...user_post\n");
#endif
    }

#endif

    /* register user functions */
    /*************************************************************************************/

    void user_register_functions(void *handle)
    {
        REGISTER_RCFD_UDF( user_set_solver_dict );
        REGISTER_RCFD_UDF( user_set_phase_dict );
        REGISTER_RCFD_UDF( user_set_rec_dict );
        REGISTER_RCFD_UDF( user_set_data_dict );
        REGISTER_RCFD_UDF( user_init_data );
        REGISTER_RCFD_UDF( user_access_data_before_phase_loop );
        REGISTER_RCFD_UDF( user_post );
        
        REGISTER_RCFD_UDF( user_set_cell_dict );
        REGISTER_RCFD_UDF( user_set_face_dict );
    }

#endif
