/******************************************************************************
 * @internal
 * @file        rCFD_prepare.h
 * @copyright   (C) 2025
 *              Particulate Flow Modelling,
 *              Johannes Kepler University, Linz, Austria
 * @website     www.particulate-flow.at
 * @authors     Stefan Pirker
 * @authors     Daniel Queteschiner
 * @authors     Hannes Lumetzberger
 *****************************************************************************/

#ifndef RCFD_PREPARE_H
#define RCFD_PREPARE_H

#include "rCFD_types.h"
#include "rCFD_macros.h"

#ifndef DEBUG_INTERPOLATION
#define DEBUG_INTERPOLATION 0
#endif

/**
 * @brief Check the error between the displacement calculated based on
 *        cell-to-cell shifts and the displacement calculated via CFD.
 *
 * @param rcfd      pointer to rCFD instance
 * @param i_phase   specifies the phase for which the displacement is tested
 * @param post      flag to indicate whether function is called pre or post
 *                  interpolation of cell-to-cell shifts
 **/
void prep_test_displacement(rCFD_type *rcfd, int i_phase, int post)
{
#if RP_NODE
    int i_dim, i_cell, i_tracer;
    int i_layer = rcfd->Solver->current_layer;
    double  total_error, total_volume;
    double  **test_displ = NULL;

    // alloc and init test_displ fields
    {
        test_displ = (double**)malloc(_CELL_DICT(i_layer).number_of_cells * sizeof(double*));

        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

            test_displ[i_cell] = (double*)malloc(3 * sizeof(double));

            for(i_dim = 0; i_dim < 3; i_dim++){

                test_displ[i_cell][i_dim] = 0.0;
            }

            _C(i_layer).weight_after_shift[i_cell] = 0.0;
        }
    }

    // loop C2Cs and set displacement field
    {
        int node0, node1, c0, c1;
        double w0;

        for(i_tracer = 0; i_tracer < rcfd->Tracer->monitor_counter[i_phase]; i_tracer++){

            node0 = rcfd->Tracer->shifts[i_phase][i_tracer].node0;
            node1 = rcfd->Tracer->shifts[i_phase][i_tracer].node1;

            if(node0 == node1){

                w0 = rcfd->Tracer->shifts[i_phase][i_tracer].w0;

                if(fabs(w0) > 0.0){

                    c0 = rcfd->Tracer->shifts[i_phase][i_tracer].c0;
                    c1 = rcfd->Tracer->shifts[i_phase][i_tracer].c1;

                    for(i_dim = 0; i_dim < 3; i_dim++){

                        test_displ[c1][i_dim] =
                            (w0 * _C(i_layer).x[c0][i_dim] + _C(i_layer).weight_after_shift[c1] * test_displ[c1][i_dim])/
                            (w0 +  _C(i_layer).weight_after_shift[c1]);
                    }

                    _C(i_layer).weight_after_shift[c1] += w0;
                }
            }
        }
    }

    // analyze global error, set new displacement field
    {
        total_volume = 0.0;
        total_error =  0.0;

        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

            if(_C(i_layer).vof[0][i_cell][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

                total_volume += _C(i_layer).volume[i_cell];

                for(i_dim = 0; i_dim < 3; i_dim++){

                    if(_C(i_layer).weight_after_shift[i_cell] > 0.0){

                        total_error += fabs(_C(i_layer).volume[i_cell]
                                            * (_C(i_layer).cfd_displ[i_phase][i_cell][i_dim] - test_displ[i_cell][i_dim])
                                            / (_C(i_layer).x_max_global[i_dim] - _C(i_layer).x_min_global[i_dim])); //error independent of bounding-box-size
                    }
                }
            }
        }

        total_volume = PRF_GRSUM1(total_volume);
        total_error =  PRF_GRSUM1(total_error);

        if(total_volume > 0.0) total_error /= total_volume;
    }

#if !defined(RCFD_STANDALONE) && DEBUG_INTERPOLATION
    // write debug information to UDMI
    {
        int i_fluent_cell, i_thread;
        Domain  *d = Get_Domain(1);
        Thread  *t;

        thread_loop_c(t,d){ begin_c_loop(i_fluent_cell, t){

            i_thread = Fluent_Topo->global_cell_thread_id[THREAD_ID(t)];
            i_cell = Fluent_Topo->global_cell_id[i_thread][i_fluent_cell];

            for(i_dim = 0; i_dim < 3; i_dim++){

                if(_C(i_layer).weight_after_shift[i_cell] > 0.0){

                    C_UDMI(i_fluent_cell, t, ((post? u_c2c1_x  : u_c2c0_x ) + i_dim))  = test_displ[i_cell][i_dim];
                    C_UDMI(i_fluent_cell, t, ((post? u_c2c1_dx : u_c2c0_dx) + i_dim))  = _C(i_layer).x[i_cell][i_dim];
                    C_UDMI(i_fluent_cell, t, ((post? u_c2c1_dx : u_c2c0_dx) + i_dim)) -= test_displ[i_cell][i_dim];
                }
                else{

                    C_UDMI(i_fluent_cell, t, ((post? u_c2c1_x  : u_c2c0_x ) + i_dim)) = 0.0;
                    C_UDMI(i_fluent_cell, t, ((post? u_c2c1_dx : u_c2c0_dx) + i_dim)) = 0.0;
                }
            }

        }end_c_loop(i_fluent_cell, t)}
    }
#endif // DEBUG_INTERPOLATION

    // free local vars
    {
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

            free(test_displ[i_cell]);
        }

        free(test_displ);
    }

    report_info0(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                 "prep_test_displacement (layer %d phase %d):"
                 " %s interpolation error (intra-process approx.): %e (m/m3)",
                 i_layer, i_phase, post?"post":"pre", total_error);
#endif // RP_NODE
}

/**
 * @brief Interpolate cross-partition cell-to-cell shifts
 *
 * @param rcfd      pointer to rCFD instance
 * @param i_layer   specifies the cell layer
 * @param i_phase   specifies the phase
 * @param shift_to_me   array to be filled with cell-to-cell shifts to this compute node
 * @return the number of cell-to-cell shifts to this compute node
 */
int prep_interpolate_parallel_shifts(rCFD_type *rcfd, int i_layer, int i_phase, C2C_shift_type **shift_to_me)
{
#if RP_NODE
    int i_dim, i_node, i_shift;
    int i_shift_new, i_shift_to_my_node, i_shift_from_my_node;
    int node0, node1, c1;
    int index_in_global_list, global_inter_node_shift_count;
    int number_of_shifts_wo_calibration, shift_to_my_node_count, shifts_from_my_node_count;
    int *list_of_int = NULL, *node_to_global_shift_count = NULL, *global_to_node_shift_count = NULL;
    double *shift_to_my_node_cal_coords = NULL, *global_inter_node_shift_cal_coords = NULL;
    double *shifts_from_my_node_cal_coords = NULL, *list_of_double = NULL, *list_of_w0 = NULL;
    C2C_shift_type *shift_to_my_node = NULL, *global_inter_node_shift = NULL, *shifts_from_my_node = NULL;

    number_of_shifts_wo_calibration = rcfd->Tracer->monitor_counter[i_phase];

    // allocate and collect shift_to_my_node, shift_to_my_node_cal_coords
    {
        // get number of per-node shift_to_my_node
        {
            i_shift_to_my_node = 0;

            for(i_shift = 0; i_shift < number_of_shifts_wo_calibration; i_shift++){

                node0 = rcfd->Tracer->shifts[i_phase][i_shift].node0;
                node1 = rcfd->Tracer->shifts[i_phase][i_shift].node1;

                if((node1 == myid) && (node0 != myid)){

                    i_shift_to_my_node++;
                }
            }

            shift_to_my_node_count = i_shift_to_my_node;
        }

        // allocate lists of shift_to_my_node and shift_to_my_node_cal_coords
        if(shift_to_my_node_count > 0){

            // list of shift_to_my_node should also host new shifts after interpolation (hence 4 x shift_to_my_node_count)
            shift_to_my_node = (C2C_shift_type*)malloc((4 * shift_to_my_node_count) * sizeof(C2C_shift_type));

            // shift_to_my_node_cal_coords is only needed for uncalibrated shifts (3 x because of 3-dim coords)
            shift_to_my_node_cal_coords = (double*)malloc((3 * shift_to_my_node_count) * sizeof(double));
        }

        // fill lists of shift_to_my_node and shift_to_my_node_cal_coords
        {
            i_shift_to_my_node = 0;

            for(i_shift = 0; i_shift < number_of_shifts_wo_calibration; i_shift++){

                node0 = rcfd->Tracer->shifts[i_phase][i_shift].node0;
                node1 = rcfd->Tracer->shifts[i_phase][i_shift].node1;

                if((node1 == myid) && (node0 != myid)){

                    shiftcpy(&shift_to_my_node[i_shift_to_my_node], &rcfd->Tracer->shifts[i_phase][i_shift]);

                    c1 = shift_to_my_node[i_shift_to_my_node].c1;

                    for(i_dim = 0; i_dim < 3; i_dim++){

                        shift_to_my_node_cal_coords[(3 * i_shift_to_my_node + i_dim)] = _C(i_layer).cfd_displ[i_phase][c1][i_dim];
                    }

                    i_shift_to_my_node++;
                }
            }
        }
    }

    // nodes send shift_to_my_node to node-0
    if(myid > 0){

        PRF_CSEND_INT(node_zero, &shift_to_my_node_count, 1, myid);

        if(shift_to_my_node_count > 0){

            list_of_int = (int*)malloc(4 * shift_to_my_node_count * sizeof(int));
            list_of_w0 = (double*)malloc(shift_to_my_node_count * sizeof(double));

            for(i_shift = 0; i_shift < shift_to_my_node_count; i_shift++){

                list_of_int[(4* i_shift + 0)] = shift_to_my_node[i_shift].c0;
                list_of_int[(4* i_shift + 1)] = shift_to_my_node[i_shift].node0;
                list_of_int[(4* i_shift + 2)] = shift_to_my_node[i_shift].c1;
                list_of_int[(4* i_shift + 3)] = shift_to_my_node[i_shift].node1;

                list_of_w0[i_shift] = shift_to_my_node[i_shift].w0;
            }

            PRF_CSEND_INT(node_zero, list_of_int, 4 * shift_to_my_node_count, myid);
            PRF_CSEND_REAL(node_zero, shift_to_my_node_cal_coords, 3 * shift_to_my_node_count, myid);
            PRF_CSEND_REAL(node_zero, list_of_w0, shift_to_my_node_count, myid);

            free(list_of_int);
            free(list_of_w0);
        }
    }

    // node-0 receives shifts and collects them into global list
    if(myid == 0){

        // get node_to_global_shift_count
        {
            global_inter_node_shift_count = shift_to_my_node_count; // node-0
            node_to_global_shift_count = (int*)malloc(compute_node_count * sizeof(int));
            node_to_global_shift_count[0] = shift_to_my_node_count;

            for(i_node = 1; i_node < compute_node_count; i_node++){

                PRF_CRECV_INT(i_node, &node_to_global_shift_count[i_node], 1, i_node);
                global_inter_node_shift_count += node_to_global_shift_count[i_node];
            }
        }

        // allocate global lists
        if(global_inter_node_shift_count > 0){

            // assume interpolation adds three new c2c shifts per existing shift (hence 4 x)
            global_inter_node_shift = (C2C_shift_type*)malloc((4 * global_inter_node_shift_count) * sizeof(C2C_shift_type));
            global_inter_node_shift_cal_coords = (double*)malloc((3 * global_inter_node_shift_count) * sizeof(double));
        }

        // fill node-0 shifts into global list
        {
            for(i_shift = 0; i_shift < node_to_global_shift_count[0]; i_shift++){

                shiftcpy(&global_inter_node_shift[i_shift], &shift_to_my_node[i_shift]);

                for(i_dim = 0; i_dim < 3; i_dim++){

                    global_inter_node_shift_cal_coords[(3 * i_shift + i_dim)] = shift_to_my_node_cal_coords[(3 * i_shift + i_dim)];
                }
            }

            index_in_global_list = node_to_global_shift_count[0];
        }

        // loop nodes, receive shifts and fill global list
        {

            for(i_node = 1; i_node < compute_node_count; i_node++){

                if(node_to_global_shift_count[i_node] <= 0) continue;

                list_of_int = (int*)malloc(4 * node_to_global_shift_count[i_node] * sizeof(int));

                PRF_CRECV_INT(i_node, list_of_int, 4 * node_to_global_shift_count[i_node], i_node);
                PRF_CRECV_REAL(i_node, &global_inter_node_shift_cal_coords[3 * index_in_global_list], 3 * node_to_global_shift_count[i_node], i_node);

                list_of_w0 = (double*)malloc(node_to_global_shift_count[i_node] * sizeof(double));

                PRF_CRECV_REAL(i_node, list_of_w0, node_to_global_shift_count[i_node], i_node);

                for(i_shift = 0; i_shift < node_to_global_shift_count[i_node]; i_shift++){

                    global_inter_node_shift[(index_in_global_list + i_shift)].c0 =     list_of_int[(4* i_shift + 0)];
                    global_inter_node_shift[(index_in_global_list + i_shift)].node0 =  list_of_int[(4* i_shift + 1)];
                    global_inter_node_shift[(index_in_global_list + i_shift)].c1 =     list_of_int[(4* i_shift + 2)];
                    global_inter_node_shift[(index_in_global_list + i_shift)].node1 =  list_of_int[(4* i_shift + 3)];

                    global_inter_node_shift[(index_in_global_list + i_shift)].w0 =     list_of_w0[i_shift];

                }

                free(list_of_int);
                free(list_of_w0);

                index_in_global_list += node_to_global_shift_count[i_node];
            }
        }
    }

    // node-0 sends interpolation requests to sending nodes
    if(myid == 0){

        // get global_to_node_shift_count list and allocate node-0 lists of shifts_from_my_node and shifts_from_my_node_cal_coords
        {
            global_to_node_shift_count = (int*)malloc(compute_node_count * sizeof(int));

            for(i_node = 0; i_node < compute_node_count; i_node++){

                global_to_node_shift_count[i_node] = 0;
            }

            for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                node0 = global_inter_node_shift[i_shift].node0;

                if((node0 < 0)||(node0 > node_last)){

                    report_error(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                                 "[%d] prep_interpolate_parallel_shifts: inter node shifts: node0 %d out of range",
                                 myid, node0);
                }
                else{

                    for(i_node = 0; i_node < compute_node_count; i_node++){

                        if(i_node == node0){

                            global_to_node_shift_count[i_node]++;
                        }
                    }
                }
            }

            shifts_from_my_node_count = global_to_node_shift_count[0];

            if(shifts_from_my_node_count > 0){

                shifts_from_my_node = (C2C_shift_type*)malloc((4 * shifts_from_my_node_count)* sizeof(C2C_shift_type));
                shifts_from_my_node_cal_coords = (double*)malloc(3 * shifts_from_my_node_count * sizeof(double));
            }
        }

        // generate new node-0 shifts
        {
            i_shift_from_my_node = 0;

            for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                if(global_inter_node_shift[i_shift].node0 == 0){

                    shiftcpy(&shifts_from_my_node[i_shift_from_my_node], &global_inter_node_shift[i_shift]);

                    for(i_dim = 0; i_dim < 3; i_dim++){

                        shifts_from_my_node_cal_coords[3 * i_shift_from_my_node + i_dim] = global_inter_node_shift_cal_coords[(3 * i_shift + i_dim)];
                    }

                    i_shift_from_my_node++;
                }
            }
        }

        // generate new node>0 shifts and send them
        {
            for(i_node = 1; i_node < compute_node_count; i_node++){

                PRF_CSEND_INT(i_node, &global_to_node_shift_count[i_node], 1, myid);

                if(global_to_node_shift_count[i_node] <= 0) continue;

                list_of_int = (int*)malloc(4 * global_to_node_shift_count[i_node] * sizeof(int));
                list_of_double = (double*)malloc(3 * global_to_node_shift_count[i_node] * sizeof(double));
                list_of_w0 = (double*)malloc(global_to_node_shift_count[i_node] * sizeof(double));

                i_shift_from_my_node = 0;

                for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                    if(global_inter_node_shift[i_shift].node0 == i_node){

                        list_of_int[(4 * i_shift_from_my_node + 0)] = global_inter_node_shift[i_shift].c0;
                        list_of_int[(4 * i_shift_from_my_node + 1)] = global_inter_node_shift[i_shift].node0;
                        list_of_int[(4 * i_shift_from_my_node + 2)] = global_inter_node_shift[i_shift].c1;
                        list_of_int[(4 * i_shift_from_my_node + 3)] = global_inter_node_shift[i_shift].node1;

                        list_of_w0[i_shift_from_my_node] = global_inter_node_shift[i_shift].w0;

                        for(i_dim = 0; i_dim < 3; i_dim++){

                            list_of_double[(3 * i_shift_from_my_node + i_dim)] = global_inter_node_shift_cal_coords[(3 * i_shift + i_dim)];
                        }

                        i_shift_from_my_node++;
                    }
                }

                PRF_CSEND_INT(i_node, list_of_int, 4 * global_to_node_shift_count[i_node], myid);
                PRF_CSEND_REAL(i_node, list_of_double, 3 * global_to_node_shift_count[i_node], myid);
                PRF_CSEND_REAL(i_node, list_of_w0, global_to_node_shift_count[i_node], myid);

                free(list_of_int);
                free(list_of_double);
                free(list_of_w0);
            }
        }
    }

    // nodes receive interpolation requests
    if(myid > 0){

        PRF_CRECV_INT(node_zero, &shifts_from_my_node_count, 1, node_zero);

        if(shifts_from_my_node_count > 0){

            shifts_from_my_node = (C2C_shift_type*)malloc((4 * shifts_from_my_node_count)* sizeof(C2C_shift_type));
            shifts_from_my_node_cal_coords = (double*)malloc(3 * shifts_from_my_node_count * sizeof(double));
            list_of_int = (int*)malloc(4 * shifts_from_my_node_count * sizeof(int));
            list_of_w0 = (double*)malloc(shifts_from_my_node_count * sizeof(double));

            PRF_CRECV_INT(node_zero, list_of_int, 4 * shifts_from_my_node_count, node_zero);
            PRF_CRECV_REAL(node_zero, shifts_from_my_node_cal_coords, 3 * shifts_from_my_node_count, node_zero);
            PRF_CRECV_REAL(node_zero, list_of_w0, shifts_from_my_node_count, node_zero);

            for(i_shift = 0; i_shift < shifts_from_my_node_count; i_shift++){

                shifts_from_my_node[i_shift].c0 =     list_of_int[(4* i_shift + 0)];
                shifts_from_my_node[i_shift].node0 =  list_of_int[(4* i_shift + 1)];
                shifts_from_my_node[i_shift].c1 =     list_of_int[(4* i_shift + 2)];
                shifts_from_my_node[i_shift].node1 =  list_of_int[(4* i_shift + 3)];
                shifts_from_my_node[i_shift].w0 =     list_of_w0[i_shift];
            }

            free(list_of_int);
            free(list_of_w0);
        }
    }

    // all nodes produce new shifts (if they receive inter-node primary shifts)
    if(shifts_from_my_node_count > 0){

        //make list 2D for function parameter
        int i;
        double **shift_coords = (double**)malloc(shifts_from_my_node_count * sizeof(double*));

        for(i=0; i<shifts_from_my_node_count; i++){

            shift_coords[i] = &shifts_from_my_node_cal_coords[3*i];
        }

        int number_of_new_shifts = shift_calibrate(rcfd,
                                                   i_layer,
                                                   i_phase,
                                                   shifts_from_my_node,
                                                   shift_coords,
                                                   shifts_from_my_node_count);

        if(number_of_new_shifts > 0){

            shifts_from_my_node_count += number_of_new_shifts;
        }

        free(shift_coords);
    }

    // nodes send new shifts to node-0
    if(myid > 0){

        PRF_CSEND_INT(node_zero, &shifts_from_my_node_count, 1, myid);

        if(shifts_from_my_node_count > 0){

            list_of_int = (int*)malloc(4 * shifts_from_my_node_count * sizeof(int));
            list_of_w0 = (double*)malloc(shifts_from_my_node_count * sizeof(double));

            for(i_shift = 0; i_shift < shifts_from_my_node_count; i_shift++){

                list_of_int[(4 * i_shift + 0)] = shifts_from_my_node[i_shift].c0;
                list_of_int[(4 * i_shift + 1)] = shifts_from_my_node[i_shift].node0;
                list_of_int[(4 * i_shift + 2)] = shifts_from_my_node[i_shift].c1;
                list_of_int[(4 * i_shift + 3)] = shifts_from_my_node[i_shift].node1;
                list_of_w0[i_shift] = shifts_from_my_node[i_shift].w0;
            }

            PRF_CSEND_INT(node_zero, list_of_int, 4 * shifts_from_my_node_count, myid);
            PRF_CSEND_REAL(node_zero, list_of_w0, shifts_from_my_node_count, myid);

            free(list_of_int);
            free(list_of_w0);
        }
    }

    // node-0 receives new shifts and collects them into new global list
    if(myid == 0){

        // node-0 fills global list
        for(i_shift = 0; i_shift < shifts_from_my_node_count; i_shift++){

            shiftcpy(&global_inter_node_shift[i_shift], &shifts_from_my_node[i_shift]);
        }

        index_in_global_list = shifts_from_my_node_count;
        node_to_global_shift_count[0] = shifts_from_my_node_count;

        // node-0 receives shifts and fills global list
        for(i_node = 1; i_node < compute_node_count; i_node++){

            PRF_CRECV_INT(i_node, &node_to_global_shift_count[i_node], 1, i_node);

            if(node_to_global_shift_count[i_node] <= 0) continue;

            list_of_int = (int*)malloc(4 * node_to_global_shift_count[i_node] * sizeof(int));
            list_of_w0 = (double*)malloc(node_to_global_shift_count[i_node] * sizeof(double));

            PRF_CRECV_INT(i_node, list_of_int, 4 * node_to_global_shift_count[i_node], i_node);
            PRF_CRECV_REAL(i_node, list_of_w0, node_to_global_shift_count[i_node], i_node);

            for(i_shift = 0; i_shift < node_to_global_shift_count[i_node]; i_shift++){

                global_inter_node_shift[(index_in_global_list + i_shift)].c0    = list_of_int[(4 * i_shift + 0)];
                global_inter_node_shift[(index_in_global_list + i_shift)].node0 = list_of_int[(4 * i_shift + 1)];
                global_inter_node_shift[(index_in_global_list + i_shift)].c1    = list_of_int[(4 * i_shift + 2)];
                global_inter_node_shift[(index_in_global_list + i_shift)].node1 = list_of_int[(4 * i_shift + 3)];
                global_inter_node_shift[(index_in_global_list + i_shift)].w0    = list_of_w0[i_shift];
            }

            index_in_global_list += node_to_global_shift_count[i_node];

            free(list_of_int);
            free(list_of_w0);
        }

        global_inter_node_shift_count = index_in_global_list;
    }

    // node-0 sends new shifts to nodes
    if(myid == 0){

        for(i_node = 0; i_node < compute_node_count; i_node++){

            global_to_node_shift_count[i_node] = 0;
        }

        for(i_node = 0; i_node < compute_node_count; i_node++){

            for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                if(global_inter_node_shift[i_shift].node1 == i_node)

                    global_to_node_shift_count[i_node]++;
            }
        }

        shift_to_my_node_count = global_to_node_shift_count[0];

        // generate new node-0 shifts
        {
            i_shift_to_my_node = 0;

            for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                if(global_inter_node_shift[i_shift].node1 == 0){

                    shiftcpy(&shift_to_my_node[i_shift_to_my_node], &global_inter_node_shift[i_shift]);

                    i_shift_to_my_node++;
                }
            }
        }

        // generate new node>0 shifts and send them
        {
            for(i_node = 1; i_node < compute_node_count; i_node++){

                PRF_CSEND_INT(i_node, &global_to_node_shift_count[i_node], 1, myid);

                if(global_to_node_shift_count[i_node] <= 0) continue;

                list_of_int = (int*)malloc(4 * global_to_node_shift_count[i_node] * sizeof(int));
                list_of_w0 = (double*)malloc(global_to_node_shift_count[i_node] * sizeof(double));

                i_shift_new = 0;

                for(i_shift = 0; i_shift < global_inter_node_shift_count; i_shift++){

                    if(global_inter_node_shift[i_shift].node1 == i_node){

                        list_of_int[(4 * i_shift_new + 0)] = global_inter_node_shift[i_shift].c0;
                        list_of_int[(4 * i_shift_new + 1)] = global_inter_node_shift[i_shift].node0;
                        list_of_int[(4 * i_shift_new + 2)] = global_inter_node_shift[i_shift].c1;
                        list_of_int[(4 * i_shift_new + 3)] = global_inter_node_shift[i_shift].node1;
                        list_of_w0[i_shift_new] = global_inter_node_shift[i_shift].w0;

                        i_shift_new++;
                    }
                }

                PRF_CSEND_INT(i_node, list_of_int, 4 * global_to_node_shift_count[i_node], myid);
                PRF_CSEND_REAL(i_node, list_of_w0, global_to_node_shift_count[i_node], myid);

                free(list_of_int);
                free(list_of_w0);
            }
        }
    }

    // nodes receive interpolation requests
    if(myid > 0){

        PRF_CRECV_INT(node_zero, &shift_to_my_node_count, 1, node_zero);

        if(shift_to_my_node_count > 0){

            list_of_int = (int*)malloc(4 * shift_to_my_node_count * sizeof(int));
            list_of_w0 = (double*)malloc(shift_to_my_node_count * sizeof(double));

            PRF_CRECV_INT(node_zero, list_of_int, 4 * shift_to_my_node_count, node_zero);
            PRF_CRECV_REAL(node_zero, list_of_w0, shift_to_my_node_count, node_zero);

            for(i_shift = 0; i_shift < shift_to_my_node_count; i_shift++){

                shift_to_my_node[i_shift].c0 =     list_of_int[(4* i_shift + 0)];
                shift_to_my_node[i_shift].node0 =  list_of_int[(4* i_shift + 1)];
                shift_to_my_node[i_shift].c1 =     list_of_int[(4* i_shift + 2)];
                shift_to_my_node[i_shift].node1 =  list_of_int[(4* i_shift + 3)];
                shift_to_my_node[i_shift].w0 =     list_of_w0[i_shift];
            }

            free(list_of_int);
            free(list_of_w0);
        }
    }

    {
        free(shift_to_my_node_cal_coords);
        free(shifts_from_my_node);
        free(shifts_from_my_node_cal_coords);

        if(myid == 0){

            free(node_to_global_shift_count);
            free(global_to_node_shift_count);
            free(global_inter_node_shift);
            free(global_inter_node_shift_cal_coords);
        }
    }

    *shift_to_me = shift_to_my_node;

    return shift_to_my_node_count;
#endif // RP_NODE
    return 0;
}

/**
 * @brief This function calculates the standard deviation from the cell-center positions
 *        of a number of shifts and returns the corresponding diffusivity of these shifts
 *
 * @param new_shifts        array of C2C shifts whose standard deviation of the origin points should be calculated
 * @param len_new_shifts    length of new_shifts array
 * @param timestep_size     used as a scaling parameter in the calculation of diffusivity
 * @param cell_centers      list of all cell_centers that are used to calculate the standard deviation.
 * @return the diffusivity of the given array of cell-to-cell shifts via their standard deviation
 */
double prep_get_interpolation_diffusivity(const C2C_shift_type *new_shifts, int len_new_shifts, double timestep_size, double **cell_centers)
{
    int i, i_dim;
    double diff = 0.0, n = 0.0;
    double avg_displ[3] = {0.0};

    // calc number of valid shifts and store in n
    for(i = 0; i < len_new_shifts; i++){

        if(new_shifts[i].w0 > 0.0){

            n++;
        }
    }

    //calc standard deviation of shift origins
    for(i_dim = 0; i_dim < 3; i_dim++){

        // calc average of shift origins
        for(i = 0; i < len_new_shifts; i++){

            if(new_shifts[i].w0 > 0.0){

                avg_displ[i_dim] += new_shifts[i].w0 * cell_centers[i][i_dim];
            }
        }

        // calc average deviation from average of shift origins
        for(i = 0; i < len_new_shifts; i++){

            if(new_shifts[i].w0 > 0.0){

                diff += new_shifts[i].w0 * (cell_centers[i][i_dim] - avg_displ[i_dim])
                                         * (cell_centers[i][i_dim] - avg_displ[i_dim]);
            }
        }
    }

    // return diffusivity
    if(n > 1)
        return (n / (n-1.0) * diff) / (2.0 * timestep_size);

    return 0;
}

/**
 * @brief count cross-partition shift cells and allocate and set outbound cell ids for each compute node
 *
 * @param rcfd          pointer to rCFD instance
 * @param requests      array obtaining number of outbound shift cells for each compute node ([i_node])
 * @param shifts        array of cross-partition cell-to-cell shifts
 * @param len_shifts    length of \p shifts array
 * @return pointer to an allocated 2d array containing outbound cell ids ([i_node][i_request])
 */
int **prep_malloc_outbound_cell_ids_per_node(rCFD_type *rcfd, int *requests, const C2C_shift_type *shifts, int len_shifts)
{
    int i_shift, j_shift, i_node, i_cell;
    int **outbound_cell_ids = (int**)malloc(compute_node_count * sizeof(int*));

    // get number of outbound requests for each node
    for(i_shift = 0; i_shift < len_shifts; i_shift++){

        if(shifts[i_shift].node0 == myid) continue;

        int is_duplicate = 0;

        for(j_shift = 0; j_shift < i_shift && !is_duplicate; j_shift++){

            if(shifts[i_shift].node0 == shifts[j_shift].node0 && shifts[i_shift].c0 == shifts[j_shift].c0){

                is_duplicate = 1;
            }
        }

        if(!is_duplicate){

            requests[shifts[i_shift].node0]++;
        }
    }

    // allocate space for storing cell_ids
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(requests[i_node] == 0){

            outbound_cell_ids[i_node] = NULL;
            continue;
        }

        outbound_cell_ids[i_node] = (int*)malloc(requests[i_node] * sizeof(int));

        for(i_cell = 0; i_cell < requests[i_node]; i_cell++){

            outbound_cell_ids[i_node][i_cell] = -1;
        }
    }

    // get requested cell_ids for each node
    for(i_shift = 0; i_shift < len_shifts; i_shift++){

        int first_free_slot = 0;

        while(outbound_cell_ids[shifts[i_shift].node0][first_free_slot] >= 0
              && outbound_cell_ids[shifts[i_shift].node0][first_free_slot] != shifts[i_shift].c0){

            first_free_slot++;
        }

        if(outbound_cell_ids[shifts[i_shift].node0][first_free_slot] < 0){ // reduntant: && first_free_slot < requests[shifts[i_shift].node0]){

            outbound_cell_ids[shifts[i_shift].node0][first_free_slot] = shifts[i_shift].c0;
        }
    }

    return outbound_cell_ids;
}

/**
 * @brief allocate space for storing inbound cell ids for each compute node according to \p requests array
 *
 * @param rcfd      pointer to rCFD instance
 * @param requests  array containing number of inbound shift cells for each compute node ([i_node])
 * @return pointer to an allocated 2d array to store inbound cell ids ([i_node][i_request])
 */
int **prep_malloc_inbound_cell_ids_per_node(rCFD_type *rcfd, const int *requests)
{
    int i_node;
    int **inbound_cell_ids = (int**)malloc(compute_node_count * sizeof(int*));

    // allocate space for storing cell_ids
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(requests[i_node] == 0){

            inbound_cell_ids[i_node] = NULL;
            continue;
        }

        inbound_cell_ids[i_node] = (int*)malloc(requests[i_node] * sizeof(int));
    }

    return inbound_cell_ids;
}

/**
 * @brief allocate space for storing a number of cell centers for each compute node as requested
 *
 * @param rcfd      pointer to rCFD instance
 * @param requests  array containing for each compute node the number of cell centers to allocate
 * @return a pointer to an allocated 3d array [i_node][i_cell][i_dim] for storing the requested number of cell centers
 */
double ***prep_malloc_cell_centers(rCFD_type *rcfd, const int *requests)
{
    int i_node;
    double ***cell_centers_per_node = (double***)malloc(compute_node_count * sizeof(double**));

    // allocate space for storing cell_centers
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(requests[i_node] == 0){

            cell_centers_per_node[i_node] = NULL;
            continue;
        }

        cell_centers_per_node[i_node] = malloc_r_2d(requests[i_node], 3);
    }

    return cell_centers_per_node;
}

/**
 * @brief allocate and set cell centers for cross-partition cell-to-cell shifts
 *
 * @param rcfd          pointer to rCFD instance
 * @param shifts        array of cross-partition cell-to-cell shifts
 * @param len_shifts    length of \p shifts array
 * @param local_cell_centers array of local cell centers
 * @return a pointer to an allocated 2d array [i_shift][i_dim] containing cell centers related to the given cell-to-cell shifts
 */
double **prep_malloc_cell_centers_parallel(rCFD_type *rcfd, const C2C_shift_type *shifts, int len_shifts, double **local_cell_centers)
{
    int i_node, i_shift, i_cell, i_dim;
#if defined(RCFD_STANDALONE)
    MPI_Request requests[2*compute_node_count];
#endif

    // set up the cell_centers to be returned to the caller
    double **cell_centers = NULL;

    if(len_shifts > 0){

        cell_centers = malloc_r_2d(len_shifts, 3);
    }

    // malloc and init requests_per_node
    int *inbound_requests_per_node = (int*)malloc(compute_node_count * sizeof(int));
    int *outbound_requests_per_node = (int*)malloc(compute_node_count * sizeof(int));

    for(i_node = 0; i_node < compute_node_count; i_node++){

        inbound_requests_per_node[i_node] = 0; //not necessary - overwritten by RECV
        outbound_requests_per_node[i_node] = 0;
    }

    // calculate the amount of outbound requests as well as the requested cell IDs for each node
    int     **outbound_cell_ids     = prep_malloc_outbound_cell_ids_per_node(rcfd, outbound_requests_per_node, shifts, len_shifts);
    double ***outbound_cell_centers = prep_malloc_cell_centers(rcfd, outbound_requests_per_node);

    // send amout of requests for cell_centers to each node
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(i_node == myid){
#if defined(RCFD_STANDALONE)
            requests[i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_ISEND_INT(i_node, &outbound_requests_per_node[i_node], 1, myid, &requests[i_node]);
#else
        PRF_ISEND_INT(i_node, &outbound_requests_per_node[i_node], 1, myid);
#endif
    }

    // receive amount of requests for cell_centers from each node
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(i_node == myid){
#if defined(RCFD_STANDALONE)
            requests[compute_node_count+i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_IRECV_INT(i_node, &inbound_requests_per_node[i_node], 1, i_node, &requests[compute_node_count+i_node]);
#else
        PRF_IRECV_INT(i_node, &inbound_requests_per_node[i_node], 1, i_node);
#endif
    }

#if defined(RCFD_STANDALONE)
    PRF_WAIT_ALL(2*compute_node_count, requests);
#else
    PRF_WAIT_ALL();
#endif

    // malloc space for the inbound cell_ids
    int     **inbound_cell_ids     = prep_malloc_inbound_cell_ids_per_node(rcfd, inbound_requests_per_node);
    double ***inbound_cell_centers = prep_malloc_cell_centers(rcfd, inbound_requests_per_node);

    // send outbound cell_ids
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(outbound_requests_per_node[i_node] == 0){
#if defined(RCFD_STANDALONE)
            requests[i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_ISEND_INT(i_node, &outbound_cell_ids[i_node][0], outbound_requests_per_node[i_node], myid, &requests[i_node]);
#else
        PRF_ISEND_INT(i_node, &outbound_cell_ids[i_node][0], outbound_requests_per_node[i_node], myid);
#endif
    }

    // receive inbound cell_ids
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(inbound_requests_per_node[i_node] == 0){
#if defined(RCFD_STANDALONE)
            requests[compute_node_count+i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_IRECV_INT(i_node, &inbound_cell_ids[i_node][0], inbound_requests_per_node[i_node], i_node, &requests[compute_node_count+i_node]);
#else
        PRF_IRECV_INT(i_node, &inbound_cell_ids[i_node][0], inbound_requests_per_node[i_node], i_node);
#endif
    }

#if defined(RCFD_STANDALONE)
    PRF_WAIT_ALL(2*compute_node_count, requests);
#else
    PRF_WAIT_ALL();
#endif

    // set up cell_centers for sending them by using the local list of cell centers
    for(i_node = 0; i_node < compute_node_count; i_node++){

        for(i_cell = 0; i_cell < inbound_requests_per_node[i_node]; i_cell++){

            for(i_dim = 0; i_dim < 3; i_dim++){

                inbound_cell_centers[i_node][i_cell][i_dim] = local_cell_centers[inbound_cell_ids[i_node][i_cell]][i_dim];
            }
        }
    }

    // send cell_centers to nodes
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(inbound_requests_per_node[i_node] == 0){
#if defined(RCFD_STANDALONE)
            requests[i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_ISEND_REAL(i_node, &inbound_cell_centers[i_node][0][0], inbound_requests_per_node[i_node] * 3, myid, &requests[i_node]);
#else
        PRF_ISEND_REAL(i_node, &inbound_cell_centers[i_node][0][0], inbound_requests_per_node[i_node] * 3, myid);
#endif
    }

    // receive cell_centers from nodes
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(outbound_requests_per_node[i_node] == 0){
#if defined(RCFD_STANDALONE)
            requests[compute_node_count+i_node] = MPI_REQUEST_NULL;
#endif
            continue;
        }

#if defined(RCFD_STANDALONE)
        PRF_IRECV_REAL(i_node, &outbound_cell_centers[i_node][0][0], outbound_requests_per_node[i_node] * 3, i_node, &requests[compute_node_count+i_node]);
#else
        PRF_IRECV_REAL(i_node, &outbound_cell_centers[i_node][0][0], outbound_requests_per_node[i_node] * 3, i_node);
#endif
    }

#if defined(RCFD_STANDALONE)
    PRF_WAIT_ALL(2*compute_node_count, requests);
#else
    PRF_WAIT_ALL();
#endif

    // store the cell centers from multiple nodes in a single list, whose entries contain the cell center of the corresponding i_shift
    for(i_node = 0; i_node < compute_node_count; i_node++){

        for(i_cell = 0; i_cell < outbound_requests_per_node[i_node]; i_cell++){

            for(i_shift = 0; i_shift < len_shifts; i_shift++){

                if(shifts[i_shift].c0 == outbound_cell_ids[i_node][i_cell] && shifts[i_shift].node0 == i_node){

                    for(i_dim = 0; i_dim < 3; i_dim++){

                        cell_centers[i_shift][i_dim] = outbound_cell_centers[i_node][i_cell][i_dim];
                    }
                }
            }
        }
    }

    // store local cell centers in the new list as well
    for(i_shift = 0; i_shift < len_shifts; i_shift++){

        if(shifts[i_shift].node0 == myid){

            for(i_dim = 0; i_dim < 3; i_dim++){

                cell_centers[i_shift][i_dim] = local_cell_centers[shifts[i_shift].c0][i_dim];
            }
        }
    }

    // free everything that was allocated except cell_centers
    for(i_node = 0; i_node < compute_node_count; i_node++){

        if(outbound_requests_per_node[i_node] > 0){

            free(outbound_cell_ids[i_node]);
            free_r_2d(outbound_cell_centers[i_node]);
        }

        if(inbound_requests_per_node[i_node] > 0){

            free(inbound_cell_ids[i_node]);
            free_r_2d(inbound_cell_centers[i_node]);
        }
    }

    free(outbound_cell_ids);
    free(inbound_cell_ids);
    free(outbound_cell_centers);
    free(inbound_cell_centers);

    free(outbound_requests_per_node);
    free(inbound_requests_per_node);

    return cell_centers;
}

/**
 * @brief set local diffusion field to diffusion stemming from interpolation
 *
 * @param local_diff    local diffusion array to be set
 * @param shifts        array of cell-to-cell shifts
 * @param len_shifts    length of \p shifts array
 * @param timestep_size global time step size
 * @param cell_centers  array of local cell centers if \p is_parallel is false, else array of external cell centers on other compute nodes related to cell-to-cell shifts
 * @param is_parallel   flag indicating if shifts are cross-partition
 */
void prep_set_local_diffusion_field(double *local_diff, const C2C_shift_type *shifts, int len_shifts, double timestep_size, double **cell_centers, int is_parallel)
{
    int i_shift_prim, i_shift, shifts_to_cell_found;
    int cells_found = 0;

    C2C_shift_type *shifts_to_cell = (C2C_shift_type*)malloc(4 * sizeof(C2C_shift_type));
    double **cell_centers_to_cell = (double**)malloc(4 * sizeof(double*));
    int *f_cell_ids;

    if(len_shifts > 0){

        f_cell_ids = (int*)malloc(len_shifts * sizeof(int));
    }
    else{

        f_cell_ids = NULL;
    }

    // make a separate list of cell_ids in order to set them to -1 when they have already been used for the diffusion calculation
    for(i_shift = 0; i_shift < len_shifts; i_shift++){

        f_cell_ids[i_shift] = shifts[i_shift].c1;
    }

    // iterate the list via i_shift_prim, find interpolation shifts via i_shift with the same ID as i_shift_prim, put them in a list, set them to -1
    i_shift_prim = 0;

    while(i_shift_prim < len_shifts){

        if(f_cell_ids[i_shift_prim] == -1){

            i_shift_prim++;
            continue;
        }

        i_shift = i_shift_prim + 1;

        shifts_to_cell[0] = shifts[i_shift_prim];

        if(is_parallel){

            cell_centers_to_cell[0] = cell_centers[i_shift_prim];
        }
        else{

            cell_centers_to_cell[0] = cell_centers[shifts[i_shift_prim].c0];
        }

        shifts_to_cell_found = 1;
        f_cell_ids[i_shift_prim] = -1;

        while(i_shift < len_shifts && shifts_to_cell_found < 4){

            if(f_cell_ids[i_shift] == shifts[i_shift_prim].c1){

                shifts_to_cell[shifts_to_cell_found] = shifts[i_shift];

                if(is_parallel){

                    cell_centers_to_cell[shifts_to_cell_found] = cell_centers[i_shift];
                }
                else{

                    cell_centers_to_cell[shifts_to_cell_found] = cell_centers[shifts[i_shift].c0];
                }

                shifts_to_cell_found++;
                f_cell_ids[i_shift] = -1;
            }

            i_shift++;
        }

        local_diff[shifts[i_shift_prim].c1] = prep_get_interpolation_diffusivity(shifts_to_cell,
                                                                                 shifts_to_cell_found,
                                                                                 timestep_size,
                                                                                 cell_centers_to_cell);

        i_shift_prim++;
        cells_found++;
    }

    free(shifts_to_cell);
    free(cell_centers_to_cell);
}

/**
 * @brief Generate additional cell-to-cell shifts via interpolation for more
 *        accurate shift patterns.
 *
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_interpolate_shifts(rCFD_type *rcfd)
{
#if RP_NODE

    /* known limitations:
        ->  single layer
    */

    // interpolation requires displacement fields
    if(rcfd->Solver_Dict->shifts_from_displacement_fields_on && rcfd->Solver_Dict->interpolation_on){

        int     i_dim, i_phase, i_shift;
        int     i_layer = rcfd->Solver->current_layer;
        int     i_per_node_shift;
        int     node0, node1, c1;
        int     number_of_shifts_wo_calibration, number_of_per_node_shifts, shift_to_my_node_count = 0;
        double  **per_node_shift_cal_coord = NULL;
        C2C_shift_type  *per_node_shift = NULL, *shift_to_my_node = NULL;

        // generate interpolation C2C shifts for all phases
        {
            for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                // propagate information to get initial displacement pattern (only intra-process)
                prep_test_displacement(rcfd, i_phase, 0);

                // create new c2c shifts by interpolation
                {
                    // initializations, allocation of per_node_shift(_cal_coords)
                    {
                        number_of_shifts_wo_calibration = rcfd->Tracer->monitor_counter[i_phase];

                        i_per_node_shift = 0;

                        for(i_shift = 0; i_shift < number_of_shifts_wo_calibration; i_shift++){

                            node0 = rcfd->Tracer->shifts[i_phase][i_shift].node0;
                            node1 = rcfd->Tracer->shifts[i_phase][i_shift].node1;

                            if(node0 == node1){

                                i_per_node_shift++;
                            }
                        }

                        number_of_per_node_shifts = i_per_node_shift;

                        if(number_of_per_node_shifts > 0){

                            // x 4 because of new interpolation shifts
                            per_node_shift = (C2C_shift_type*)malloc(4 * number_of_per_node_shifts * sizeof(C2C_shift_type));

                            // x 3 because of 3-dim coords
                            per_node_shift_cal_coord = malloc_r_2d(number_of_per_node_shifts, 3);
                        }
                    }

                    // consider per_node tracers (intra-process)
                    {
                        // fill per_node_shifts(_cal_coords)
                        {
                            i_per_node_shift = 0;

                            for(i_shift = 0; i_shift < number_of_shifts_wo_calibration; i_shift++){

                                node0 = rcfd->Tracer->shifts[i_phase][i_shift].node0;
                                node1 = rcfd->Tracer->shifts[i_phase][i_shift].node1;

                                if(node0 == node1){

                                    shiftcpy(&per_node_shift[i_per_node_shift], &rcfd->Tracer->shifts[i_phase][i_shift]);

                                    c1 = per_node_shift[i_per_node_shift].c1;

                                    for(i_dim = 0; i_dim < 3; i_dim++){

                                        per_node_shift_cal_coord[i_per_node_shift][i_dim] = _C(i_layer).cfd_displ[i_phase][c1][i_dim];    // real-world coord.
                                    }

                                    i_per_node_shift++;
                                }
                            }
                        }

                        // add new tracers by interpolation
                        if(number_of_per_node_shifts > 0){

                            int number_of_new_shifts = shift_calibrate(rcfd,
                                                                       i_layer,
                                                                       i_phase,
                                                                       per_node_shift,
                                                                       per_node_shift_cal_coord,
                                                                       number_of_per_node_shifts);

                            if(number_of_new_shifts > 0){

                                number_of_per_node_shifts += number_of_new_shifts;
                            }
                        }
                    }

                    // consider inter-process shifts
                    if(rcfd->Calibration_Dict->calibration_parallel_on){

                         shift_to_my_node_count = prep_interpolate_parallel_shifts(rcfd, i_layer, i_phase, &shift_to_my_node);
                    }

                    // refill rcfd->Tracer->shifts[i_phase] with per_node_shift and shift_to_my_node
                    {
                        for(i_shift = 0; i_shift < number_of_per_node_shifts; i_shift++){

                            shiftcpy(&rcfd->Tracer->shifts[i_phase][i_shift], &per_node_shift[i_shift]);
                        }

                        if(rcfd->Calibration_Dict->calibration_parallel_on){

                            for(i_shift = 0; i_shift < shift_to_my_node_count; i_shift++){

                                shiftcpy(&rcfd->Tracer->shifts[i_phase][(number_of_per_node_shifts + i_shift)], &shift_to_my_node[i_shift]);
                            }
                        }

                        rcfd->Tracer->monitor_counter[i_phase] = number_of_per_node_shifts + shift_to_my_node_count;
                    }

                    // free local vars
                    {
                        free(per_node_shift);
                        per_node_shift = NULL;

                        free_r_2d(per_node_shift_cal_coord);
                        per_node_shift_cal_coord = NULL;

                        free(shift_to_my_node);
                        shift_to_my_node = NULL;
                    }
                }

                // propagate information to get final displacement pattern (only intra-process)
                prep_test_displacement(rcfd, i_phase, 1);

                if(rcfd->Solver_Dict->local_diff_on){

                    // clear diffusion array
                    memset(&_C(i_layer).local_diff[0][0][0], 0, sizeof(double) * 1 * rcfd->Solver_Dict->number_of_phases * _CELL_DICT(i_layer).number_of_cells);

                    prep_set_local_diffusion_field(_C(i_layer).local_diff[0][i_phase],
                                                    rcfd->Tracer->shifts[i_phase],
                                                    number_of_per_node_shifts,
                                                    rcfd->Solver_Dict->global_time_step,
                                                    _C(i_layer).x,
                                                    0);

                    double **cell_centers_parallel = prep_malloc_cell_centers_parallel(rcfd,
                                                                                        &(rcfd->Tracer->shifts[i_phase][number_of_per_node_shifts]),
                                                                                        shift_to_my_node_count,
                                                                                        _C(i_layer).x);

                    // no need to set a diffusion field or free cell_centers_parallel if there are no parallel shifts
                    // (but the allocation process above uses parallel communication, so it always needs to happen)
                    if(shift_to_my_node_count > 0){

                        prep_set_local_diffusion_field(_C(i_layer).local_diff[0][i_phase],
                                                        &(rcfd->Tracer->shifts[i_phase][number_of_per_node_shifts]),
                                                        shift_to_my_node_count,
                                                        rcfd->Solver_Dict->global_time_step,
                                                        cell_centers_parallel,
                                                        1);
                    }

                    free_r_2d(cell_centers_parallel);
                    cell_centers_parallel = NULL;
                }
            }
        }
    }
    else{

        if (myid == 0){

            report_error(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                         "prep_interpolate_shifts: interpolation %s",
                         (rcfd->Solver_Dict->interpolation_on) ? "requires displacement fields" : "is turned off");
        }

        return -1;
    }
#endif // RP_NODE

    return 0;
}

/**
 * @brief write the recurrence matrix, recurrence path and recurrence path
 *        differences to files
 *
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_recurrence(rCFD_type *rcfd)
{
#if RP_NODE
    double  **Norm_Database;
    int     i_state, i_norm, i_cell, i_tmp, i_layer = 0;
    FILE    *f_in = NULL;
#endif // RP_NODE

    int     i_frame = 0, i_frame2 = 0, i_island;
    int     number_of_jumps;
    double  norm, diff;
    char    file_name[80];
    int success = 0;

#if RP_HOST
    FILE    *f_out = NULL;
    double  **Rec_matrix = NULL;
    int     *Rec_jump = NULL;
    double  *Rec_jump_diff = NULL;

#ifdef RCFD_STANDALONE
    if(myid == 0)
#endif
    {
        Rec_matrix      = malloc_r_2d(rcfd->Solver_Dict->number_of_frames, rcfd->Solver_Dict->number_of_frames);
        Rec_jump        = malloc(sizeof(int) * rcfd->Solver_Dict->number_of_frames);
        Rec_jump_diff   = malloc(sizeof(double) * rcfd->Solver_Dict->number_of_frames);
    }
#endif // RP_HOST

    /* read Norms */
    {
#if RP_NODE

        Norm_Database = malloc_r_2d(rcfd->Solver_Dict->number_of_frames, rcfd->Norms->number_of_norms);

        i_state = rcfd->Solver->current_state;

        sprintf(file_name,"%s_%d_%d", rcfd->File_Dict->Norm_filename, i_state, myid);

        f_in = fopen(file_name,"r");

        if(f_in == NULL){

            report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE,
                         "[%d] prep_write_recurrence: failed to open: '%s'.",
                         myid, file_name);
            success = -1;
        }
        else{

            for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                fscanf(f_in, "%d \n", &i_tmp);

                for(i_norm = 0; i_norm < rcfd->Norms->number_of_norms; i_norm++){

                    fscanf(f_in, "%le \n", &Norm_Database[i_frame][i_norm]);
                }
            }

            fclose(f_in);
        }

        success = PRF_GILOW1(success);
#endif // RP_NODE
    }

    node_to_host_int_1(success);
    if(success < 0)
        return success;

    number_of_jumps = 0;

    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

        /* calc Rec_Matrix */
        {
#if RP_NODE
            double  mean_norm, sum_of_all_norms = 0.0;

            for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                for(i_norm = 0; i_norm < rcfd->Norms->number_of_norms; i_norm++){

                    sum_of_all_norms += Norm_Database[i_frame][i_norm];
                }
            }

            mean_norm = PRF_GRSUM1(sum_of_all_norms) /
                (double) PRF_GISUM1(rcfd->Norms->number_of_norms) /
                (double) rcfd->Solver_Dict->number_of_frames;
#endif // RP_NODE

            for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                for(i_frame2 = 0; i_frame2 < rcfd->Solver_Dict->number_of_frames; i_frame2++){
#if RP_NODE
                    norm = 0.0;

                    for(i_norm = 0; i_norm < rcfd->Norms->number_of_norms; i_norm++){

                        i_cell = i_norm * rcfd->Norm_Dict->coarse_graining;

                        if( _C(i_layer).island_id[i_cell] == i_island){

                            diff = fabs(Norm_Database[i_frame][i_norm] - Norm_Database[i_frame2][i_norm]);

                            if(rcfd->Norm_Dict->larger_than_mean_norm_only){

                                if(diff > mean_norm){

                                    norm += diff;
                                }
                            }
                            else{

                                norm += diff;
                            }
                        }
                    }

                    norm = PRF_GRSUM1(norm);
#endif // RP_NODE

                    node_to_host_real_1(norm);

#if RP_HOST
#ifdef RCFD_STANDALONE
                    if(myid == 0)
#endif
                    {
                        Rec_matrix[i_frame][i_frame2] = norm;
                        Rec_matrix[i_frame2][i_frame] = norm;
                    }
#endif // RP_HOST
                }
            }
        }

        /* calc jumps */
#ifdef RCFD_STANDALONE
        if(myid == 0)
#endif
        {
#if RP_HOST
            switch (rcfd->Rec_Dict->format){

                case RCFD_REC_QUARTER_JUMPS:
                {
                    int N1_4=(int)floor((double)rcfd->Solver_Dict->number_of_frames/4.);
                    int N2_4=(int)floor((double)rcfd->Solver_Dict->number_of_frames/2.);
                    int N3_4=(int)floor((double)rcfd->Solver_Dict->number_of_frames*3./4.);

                    int j_frame, j_frame_min = 0;

                    for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                        diff = 1e10; /* something large */

                        if(i_frame <= N1_4){

                            for(j_frame = N2_4; j_frame < N3_4; j_frame++){

                                if(Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }
                        else if(i_frame <= N2_4){

                            for(j_frame = N3_4; j_frame < rcfd->Solver_Dict->number_of_frames; j_frame++){

                                if (Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }
                        else if(i_frame <= N3_4){

                            for(j_frame = 0; j_frame < N1_4; j_frame++){

                                if (Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }
                        else{

                            for(j_frame = N1_4; j_frame < N2_4; j_frame++){

                                if (Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }

                        Rec_jump[i_frame] = j_frame_min;
                        Rec_jump_diff[i_frame] = diff;
                    }

                    break;
                }
                case RCFD_REC_OFF_DIAGONAL_BAND:
                {
                    int j_frame, j_frame_min = 0;

                    for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                        diff = 1e10;

                        if(i_frame > rcfd->Rec_Dict->off_diagonal_band_width){

                            for(j_frame = 0; j_frame < (i_frame - rcfd->Rec_Dict->off_diagonal_band_width); j_frame++){

                                if (Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }

                        if(i_frame < (rcfd->Solver_Dict->number_of_frames - rcfd->Rec_Dict->off_diagonal_band_width)){

                            for(j_frame = (i_frame + rcfd->Rec_Dict->off_diagonal_band_width); j_frame < rcfd->Solver_Dict->number_of_frames; j_frame++){

                                if (Rec_matrix[i_frame][j_frame] < diff){

                                    diff = Rec_matrix[i_frame][j_frame];
                                    j_frame_min = j_frame;
                                }
                            }
                        }

                        Rec_jump[i_frame] = j_frame_min;
                        Rec_jump_diff[i_frame] = diff;
                    }

                    break;
                }
                case RCFD_REC_REPLAY:
                {
                    for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                        if(i_frame < (rcfd->Solver_Dict->number_of_frames - 1)){

                            Rec_jump[i_frame] = i_frame + 1;
                            Rec_jump_diff[i_frame] = Rec_matrix[i_frame][(i_frame + 1)];
                        }
                        else{

                            Rec_jump[i_frame] = 0;
                            Rec_jump_diff[i_frame] = Rec_matrix[i_frame][0];
                        }
                    }
                    break;
                }
            }
#endif // RP_HOST
        }

        /* write matrix, jumps */
#ifdef RCFD_STANDALONE
        if(myid == 0)
#endif
        {
#if RP_HOST
            /* Rec Matrix */
            {
                sprintf(file_name, "%s_%d_%d_%d", rcfd->File_Dict->Matrix_filename, rcfd->Solver->current_state, rcfd->Solver->current_state, i_island);

                f_out = fopen(file_name, "w");

                if(f_out == NULL){

                    report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE,
                                 "prep_write_recurrence: failed to open: '%s'.", file_name);
                    success = -1;
                    break;
                }

                for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                    for(i_frame2 = 0; i_frame2 < rcfd->Solver_Dict->number_of_frames; i_frame2++){

                        fprintf(f_out, "%e ", Rec_matrix[i_frame][i_frame2]);
                    }

                    fprintf(f_out, "\n");
                }

                report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE,
                             "prep_write_recurrence: wrote %dx%d recurrence matrix to '%s'",
                             i_frame, i_frame2, file_name);
                fclose(f_out);
            }

            /* Rec Path */
            {
                sprintf(file_name, "%s_%d_%d_%d", rcfd->File_Dict->Jump_filename, rcfd->Solver->current_state, rcfd->Solver->current_state, i_island);

                f_out = fopen(file_name, "w");

                if(f_out == NULL){

                    report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE, "prep_write_recurrence: failed to open: '%s'.", file_name);
                    success = -1;
                    break;
                }

                for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                    fprintf(f_out, "%d \n", Rec_jump[i_frame]);
                }

                report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE, "prep_write_recurrence: wrote %d path entries to '%s'", i_frame, file_name);
                fclose(f_out);

                ++number_of_jumps;
            }

            /* Rec Path Diff */
            {
                sprintf(file_name, "%s_diff_%d_%d_%d", rcfd->File_Dict->Jump_filename, rcfd->Solver->current_state, rcfd->Solver->current_state, i_island);

                f_out = fopen(file_name, "w");

                if(f_out == NULL){

                    report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE, "prep_write_recurrence: failed to open: '%s'.", file_name);
                    success = -1;
                    break;
                }

                for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                    fprintf(f_out, "%e \n", Rec_jump_diff[i_frame]);
                }

                report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE, "prep_write_recurrence: wrote %d path diff entries to '%s'", i_frame, file_name);
                fclose(f_out);
            }
#endif // RP_HOST
        }
    }

    host_to_node_int_1(success);

    /* free local vars */
    {
#if RP_NODE
        free_r_2d(Norm_Database);
#endif // RP_NODE

#if RP_HOST
        free_r_2d(Rec_matrix);
        free(Rec_jump);
        free(Rec_jump_diff);
#endif
    }

    report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_RECURRENCE, "prep_write_recurrence: wrote %d jump vector(s)", number_of_jumps);

    return success;
}

/**
 * @brief set the norms, vof, local diffusion and user-defined recurrence fields for the current frame
 * @param rcfd  pointer to rCFD instance
 * @param norm  array of cell norms; may be NULL if set directly or via user defined function
 * @param vof   vof field for each phase ([i_cell][i_phase]); may be NULL in case of single-phase or if set directly
 * @param local_diffusion  local diffusion field ([i_phase][i_cell]); may be NULL if local diffusion is turned off or if set directly
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_set_fields(rCFD_type *rcfd, double *norm, double **vof, double **local_diffusion)
{
#if RP_NODE
    int i_cell, i_phase, i_layer = 0;

    if(rcfd->Norm_Dict->format == RCFD_NORM_DEFAULT){

        if(norm){

            memcpy(rcfd->Norms->norm, norm, sizeof(double) * rcfd->Norms->number_of_norms);
        }
    }
    else{

        rcfd->rCFD_UDF->_user_set_norm(rcfd);
    }

    if((rcfd->Solver_Dict->number_of_phases > 1) || (rcfd->Solver_Dict->single_phase_vof)){

        if(vof){

            // must not assume that input array is contiguous memory so loop instead of memcpy
            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; ++i_cell){

                for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; ++i_phase){

                    _C(i_layer).vof[0][i_cell][i_phase] = vof[i_cell][i_phase];
                }
            }
        }
    }

    if(rcfd->Solver_Dict->local_diff_on){

        if(local_diffusion){

            // must not assume that input array is contiguous memory so loop instead of memcpy
            for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; ++i_phase){

                for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; ++i_cell){

                    _C(i_layer).local_diff[0][i_phase][i_cell] = local_diffusion[i_phase][i_cell];
                }
            }
        }
    }

    rcfd->rCFD_UDF->_user_set_rec_vars(rcfd, 0);

#endif // RP_NODE

    return 0;
}

/**
 * @brief write a database entry for the norms of the current frame
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_norm(rCFD_type *rcfd)
{
#if RP_NODE
    int i_norm;
    int i_state = rcfd->Solver->current_state;
    int total_number_of_norms;
    char file_name[80];
    FILE *f_out = NULL;

    sprintf(file_name,"%s_%d_%d", rcfd->File_Dict->Norm_filename, i_state, myid);

    f_out = fopen(file_name, (rcfd->Norms->frame_counter == 0) ? "w" : "a");

    if(f_out == NULL){

        report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS, "[%d] prep_write_norm: failed to open file %s", myid, file_name);
        return -1;
    }

    fprintf(f_out, "%d \n", rcfd->Norms->number_of_norms);

    for(i_norm = 0; i_norm < rcfd->Norms->number_of_norms; i_norm++){

        fprintf(f_out,"%e \n", rcfd->Norms->norm[i_norm]);
    }

    fclose(f_out);

    total_number_of_norms = PRF_GISUM1(rcfd->Norms->number_of_norms);

    report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS, "prep_write_norm: wrote norm frame # %d with %d values", rcfd->Norms->frame_counter, total_number_of_norms);

#endif // RP_NODE

    return 0;
}

/**
 * @brief write a database entry for the vof field(s) of the current frame
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_vof(rCFD_type *rcfd)
{
#if RP_NODE
    if((rcfd->Solver_Dict->number_of_phases > 1) || (rcfd->Solver_Dict->single_phase_vof)){

        int i_phase, i_cell;
        int i_layer = 0;
        int total_number_of_cells;
        int i_state = rcfd->Solver->current_state;
        char file_name[80];
        FILE *f_out = NULL;

        for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

            sprintf(file_name, "%s_%d_%d_%d", rcfd->File_Dict->vof_filename, i_state, i_phase, myid);

            f_out = fopen(file_name, (rcfd->Norms->frame_counter == 0) ? "w" : "a");

            if(f_out == NULL){

                report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                             "[%d] prep_write_vof: failed to open file %s", myid, file_name);
                return -1;
            }

            fprintf(f_out,"%d \n", _CELL_DICT(i_layer).number_of_cells);

            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                fprintf(f_out,"%e \n", _C(i_layer).vof[0][i_cell][i_phase]);
            }

            fclose(f_out);

            total_number_of_cells = PRF_GISUM1(_CELL_DICT(i_layer).number_of_cells);

            report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                         "prep_write_vof: wrote vof frame # %d of phase %d with %d values (incl. corona cells)",
                         rcfd->Norms->frame_counter, i_phase, total_number_of_cells);
        }
    }
#endif // RP_NODE

    return 0;
}

/**
 * @brief write a database entry for the local diffusion field of the current frame
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_local_diffusion(rCFD_type *rcfd)
{
#if RP_NODE
    if(rcfd->Solver_Dict->local_diff_on){

        int i_phase, i_cell;
        int total_number_of_cells;
        int i_layer = 0;
        int i_state = rcfd->Solver->current_state;
        char file_name[80];
        FILE *f_out = NULL;

        for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

            sprintf(file_name,"%s_%d_%d_%d", rcfd->File_Dict->local_diff_filename, i_state, i_phase, myid);

            f_out = fopen(file_name, (rcfd->Norms->frame_counter == 0) ? "w" : "a");

            if(f_out == NULL){

                report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                             "[%d] prep_write_local_diffusion: failed to open file %s", myid, file_name);
                return -1;
            }

            fprintf(f_out, "%d \n", _CELL_DICT(i_layer).number_of_cells);

            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    fprintf(f_out, "%e \n", _C(i_layer).local_diff[0][i_phase][i_cell]);
            }

            fclose(f_out);

            total_number_of_cells = PRF_GISUM1(_CELL_DICT(i_layer).number_of_cells);

            report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                         "prep_write_local_diffusion: wrote local_diff frame # %d of phase %d with %d values (incl. corona cells)",
                         rcfd->Norms->frame_counter, i_phase, total_number_of_cells);
        }
    }
#endif // RP_NODE

    return 0;
}

/**
 * @brief write a database entry for the user-defined recurrence field(s) of the current frame
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_rec_user(rCFD_type *rcfd)
{
#if RP_NODE
    int i_layer = 0;

    if(_CELL_DICT(i_layer).number_of_rec_user_vars > 0){

        int i_rec_user, i_cell;
        int i_state = rcfd->Solver->current_state;
        char file_name[80];
        FILE *f_out = NULL;

        sprintf(file_name, "%s_%d_%d", rcfd->File_Dict->rec_user_filename, i_state, myid);

        f_out = fopen(file_name, (rcfd->Norms->frame_counter == 0) ? "w" : "a");

        if(f_out == NULL){

            report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                         "[%d] prep_write_rec_user: failed to open file %s", myid, file_name);
            return -1;
        }

        fprintf(f_out, "%d \n", _CELL_DICT(i_layer).number_of_int_cells);

        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

            for(i_rec_user = 0; i_rec_user < _CELL_DICT(i_layer).number_of_rec_user_vars; i_rec_user++){

                 fprintf(f_out, "%e ", _C(i_layer).rec_user[0][i_cell][i_rec_user]);
             }

            fprintf(f_out, "\n");
        }

        fclose(f_out);

        report_info0(rcfd, RCFD_REPORT_CATEGORY_WRITE_FIELDS,
                     "prep_write_rec_user: frame # %d written", rcfd->Norms->frame_counter);
    }
#endif // RP_NODE

    return 0;
}

/**
 * @brief write a database entry for norms, vof, local diffusion and user-defined recurrence fields of the current frame
 * @param rcfd  pointer to rCFD instance
 * @retval -1 Error
 * @retval 0 OK
 */
int prep_write_fields(rCFD_type *rcfd)
{
    int success = 0;

#if RP_NODE

    success |= prep_write_norm(rcfd);
    success |= prep_write_vof(rcfd);
    success |= prep_write_local_diffusion(rcfd);
    success |= prep_write_rec_user(rcfd);

    rcfd->Norms->frame_counter++;

#endif // RP_NODE

    return success;
}


#endif // RCFD_PREPARE_H
