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

#ifndef RCFD_FLUENT_H
#define RCFD_FLUENT_H

#include "udf.h"
#include "rCFD_library.h"
#include "rCFD_globals.h"
#include "rCFD_types.h"

#ifndef REGISTER_RCFD_UDF
#define REGISTER_RCFD_UDF(funcName) \
        rcfd_register_user_function(handle, #funcName, (void(*)(void))funcName)
#endif

// Fluent Particle Vars
#define     p_just_started                  (0 + p_vars_used_by_case)
#define     p_start_time                    (1 + p_vars_used_by_case)
#define     p_invalid                       (2 + p_vars_used_by_case)
#define     p_phase_id                      (3 + p_vars_used_by_case)
#define     p_phase_fraction                (4 + p_vars_used_by_case)
#define     p_c0                            (5 + p_vars_used_by_case)
#define     p_node0                         (6 + p_vars_used_by_case)
#define     p_w0                            (7 + p_vars_used_by_case)
#define     p_c_old                         (8 + p_vars_used_by_case)
#define     p_time_ratio                    (9 + p_vars_used_by_case)
#define     p_u_rwm                         (10 + p_vars_used_by_case)
#define     p_v_rwm                         (11 + p_vars_used_by_case)
#define     p_w_rwm                         (12 + p_vars_used_by_case)
#define     p_vel_rwm_old                   (13 + p_vars_used_by_case)

    enum {
        NODE_VSITED = 0,
        NODE_LOCAL_ID
    };

    typedef struct  Fluent_Topo_struct
    {
        int     number_of_int_cells, number_of_ext_cells;

        int     number_of_int_faces, number_of_ext_faces;

        int     number_of_nodes;

        int     *global_cell_thread_id;

        int     **global_cell_id;

    } Fluent_Topo_type;


    /**
     * @brief allocate and initialize the global Fluent_Topo structure with the
     *        ANSYS Fluent mesh information
     */
    void fluent_init_topo()
    {
#if RP_NODE

#if 1   // local vars

        int     i_thread, i_fluent_cell, i_cell_int, i_cell_ext, i_node;

        int     max_fluent_cell_threads_id, number_of_fluent_cell_threads, max_fluent_cell_id;

        Domain  *d = Get_Domain(1);

        Thread  *t = NULL, *ct0 = NULL, *ct1 = NULL;

        Node    *node;

#endif
        Fluent_Topo = (Fluent_Topo_type*)malloc(sizeof(Fluent_Topo_type));

        // get max_fluent_cell_threads_id, set Fluent_Topo->global_cell_thread_id
        {
            max_fluent_cell_threads_id = -1;

            number_of_fluent_cell_threads = 0;

            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

                if(max_fluent_cell_threads_id < THREAD_ID(t)){

                    max_fluent_cell_threads_id = THREAD_ID(t);
                }

                number_of_fluent_cell_threads++;

            }}

            Fluent_Topo->global_cell_thread_id = (int*)malloc((max_fluent_cell_threads_id + 1) * sizeof(int));

            i_thread = 0;

            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

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

                i_thread++;

            }}
        }

        // get number of int. cells and max. of i_fluent_cell, set Fluent_Topo->global_cell_id
        {
            Fluent_Topo->number_of_int_cells = 0;

            max_fluent_cell_id = -1;

            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

                begin_c_loop_all(i_fluent_cell, t){

                    if(C_PART(i_fluent_cell, t) == myid){

                        Fluent_Topo->number_of_int_cells++;
                    }

                    if(max_fluent_cell_id < i_fluent_cell){

                        max_fluent_cell_id = i_fluent_cell;
                    }

                    c_node_loop(i_fluent_cell, t, i_node){

                        node = C_NODE(i_fluent_cell, t, i_node);

                        N_UDMI(node, NODE_VSITED) = 0.0; // not yet visited
                    }

                }end_c_loop_all(i_fluent_cell, t)

            }}

            {
#if 0
                Fluent_Topo->global_cell_id = malloc_i_2d(number_of_fluent_cell_threads, (max_fluent_cell_id + 1));
#else
                int n = 0;

                int *data = (int*)malloc(sizeof(int) * number_of_fluent_cell_threads * (max_fluent_cell_id + 1));

                Fluent_Topo->global_cell_id = (int**)malloc(sizeof(int*) * number_of_fluent_cell_threads);

                for(i_thread = 0; i_thread < number_of_fluent_cell_threads; ++i_thread){

                    Fluent_Topo->global_cell_id[i_thread] = &data[n];

                    n += (max_fluent_cell_id + 1);
                }
#endif
            }

            i_cell_int = 0;

            i_cell_ext = 0;

            Fluent_Topo->number_of_nodes = 0;

            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

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

                begin_c_loop_all(i_fluent_cell, t){

                    if(C_PART(i_fluent_cell, t) == myid){

                        Fluent_Topo->global_cell_id[i_thread][i_fluent_cell] = i_cell_int;

                        i_cell_int++;
                    }
                    else{

                        Fluent_Topo->global_cell_id[i_thread][i_fluent_cell] = Fluent_Topo->number_of_int_cells + i_cell_ext;

                        i_cell_ext++;
                    }

                    c_node_loop(i_fluent_cell, t, i_node){

                        node = C_NODE(i_fluent_cell, t, i_node);

                        if(N_UDMI(node, NODE_VSITED) == 0.0){

                            N_UDMI(node, NODE_VSITED) = 1.0; // visited

                            Fluent_Topo->number_of_nodes++;
                        }
                    }

                }end_c_loop_all(i_fluent_cell, t)

            }}

            Fluent_Topo->number_of_ext_cells = i_cell_ext;
        }

        // get number of int. faces and ext. faces
        {
            int c0, c1, i_fluent_face;

            Fluent_Topo->number_of_int_faces = 0;
            Fluent_Topo->number_of_ext_faces = 0;

            thread_loop_f(t,d){if(THREAD_TYPE(t)==THREAD_F_INTERIOR){begin_f_loop(i_fluent_face, t){

                c0 = (int)F_C0(i_fluent_face, t);

                c1 = (int)F_C1(i_fluent_face, t);

                ct0 = F_C0_THREAD(i_fluent_face, t);

                ct1 = F_C1_THREAD(i_fluent_face, t);

                if((c0 >= 0)&&(c1 >= 0)){

                    if((C_PART(c0, ct0) == myid) && (C_PART(c1, ct1) == myid)){

                        Fluent_Topo->number_of_int_faces++;
                    }
                    else if((C_PART(c0, ct0) == myid) || (C_PART(c1, ct1) == myid)){

                        Fluent_Topo->number_of_ext_faces++;
                    }
                }
            }end_f_loop(i_fluent_face, t)}}
        }
#endif // RP_NODE
    }

    /**
     * @brief free the global Fluent_Topo structure
     */
    void fluent_free_topo()
    {
#if RP_NODE
        if(Fluent_Topo){

            free(Fluent_Topo->global_cell_thread_id);

            if(Fluent_Topo->global_cell_id){

                free(Fluent_Topo->global_cell_id[0]);
                free(Fluent_Topo->global_cell_id);
            }

            free(Fluent_Topo);
        }
#endif
    }

    /**
     * @brief extract the number of internal and external cells from the ANSYS
     *        Fluent mesh and set those parameters in the rCFD Cell_Dict. Note
     *        that this overwrites the parameters that may have been read from
     *        rCFD topo files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_set_cell_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        if(i_layer == 0){

            int     number_of_int_cells, number_of_ext_cells;

            int     i_fluent_cell;

            Domain  *d = Get_Domain(1);

            Thread  *t = NULL;

            number_of_int_cells = 0;

            number_of_ext_cells = 0;

            thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){begin_c_loop_all(i_fluent_cell, t){

                if(C_PART(i_fluent_cell, t) == myid){

                    number_of_int_cells++;
                }
                else{

                    number_of_ext_cells++;
                }

            }end_c_loop_all(i_fluent_cell, t)}}

            rcfd_set_dict_property_i(handle, "Cell_Dict", i_layer, -1, -1, "number_of_int_cells", number_of_int_cells);
            rcfd_set_dict_property_i(handle, "Cell_Dict", i_layer, -1, -1, "number_of_ext_cells", number_of_ext_cells);
        }
#endif // RP_NODE
    }

    /**
     * @brief extract the number of internal and external faces from the ANSYS
     *        Fluent mesh and set those parameters in the rCFD Face_Dict. Note
     *        that this overwrites the parameters that may have been read from
     *        rCFD topo files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_set_face_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        if(i_layer == 0){

            int     number_of_int_faces, number_of_ext_faces;

            int     i_fluent_face;

            int     c0, c1;

            Domain  *d = Get_Domain(1);

            Thread  *t = NULL, *ct0 = NULL, *ct1 = NULL;

            {
                number_of_int_faces = 0;

                number_of_ext_faces = 0;

                thread_loop_f(t,d){if(THREAD_TYPE(t)==THREAD_F_INTERIOR){begin_f_loop(i_fluent_face, t){

                    c0 = (int)F_C0(i_fluent_face, t);

                    c1 = (int)F_C1(i_fluent_face, t);

                    ct0 = F_C0_THREAD(i_fluent_face, t);

                    ct1 = F_C1_THREAD(i_fluent_face, t);

                    if((c0 >= 0)&&(c1 >= 0)){

                        if((C_PART(c0, ct0) == myid) && (C_PART(c1, ct1) == myid)){

                            number_of_int_faces++;
                        }
                        else if((C_PART(c0, ct0) == myid) || (C_PART(c1, ct1) == myid)){

                            number_of_ext_faces++;
                        }
                    }

                }end_f_loop(i_fluent_face, t)}}

                rcfd_set_dict_property_i(handle, "Face_Dict", i_layer, -1, -1, "number_of_int_faces", number_of_int_faces);
                rcfd_set_dict_property_i(handle, "Face_Dict", i_layer, -1, -1, "number_of_ext_faces", number_of_ext_faces);
            }
        }
#endif // RP_NODE
    }

    /**
     * @brief set the number of cell nodes extracted from from the ANSYS Fluent
     *        mesh in the rCFD Node_Dict. Note that this overwrites the
     *        parameter that may have been read from rCFD topo files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_set_node_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        if(i_layer == 0){

            rcfd_set_dict_property_i(handle, "Node_Dict", i_layer, -1, -1, "number_of_nodes", Fluent_Topo->number_of_nodes);
        }
#endif // RP_NODE
    }

    /**
     * @brief extract cell centers and cell volumes from the ANSYS Fluent mesh
     *        to initialize the corresponding rCFD cell data structures. Note
     *        that this overwrites the values that may have been read from rCFD
     *        topo files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_init_cells(void *handle, const short i_layer)
    {
#if RP_NODE

        if(i_layer == 0){

#if 1   // local vars

            int     i_thread, i_cell, i_fluent_cell, i_dim;

            double  x[3];

            double  **cell_x    = rcfd_get_cell_center(handle, i_layer);

            double  *cell_volume= rcfd_get_cell_volume(handle, i_layer);

            Domain  *d = Get_Domain(1);

            Thread  *t;
#endif

            // set Fluent based cell values (x, vol, dx)
            {
                thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

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

                    begin_c_loop_all(i_fluent_cell, t){

                        // set _C(i_layer).x, _C(i_layer).volume
                        {
                            i_cell = Fluent_Topo->global_cell_id[i_thread][i_fluent_cell];

                            C_CENTROID(x, i_fluent_cell, t);

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

                                cell_x[i_cell][i_dim] = x[i_dim];
                            }

                            cell_volume[i_cell] = C_VOLUME(i_fluent_cell, t);
                        }

                    }end_c_loop_all(i_fluent_cell,t)
                }}
            }
        }
#endif // RP_NODE
    }

    /**
     * @brief extract face area vectors and the cells on either side of the
     *        faces from the ANSYS Fluent mesh to initialize the corresponding
     *        rCFD face data structures. Note that this overwrites the values
     *        that may have been read from rCFD topo files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_init_faces(void *handle, const short i_layer)
    {
#if RP_NODE
        if(i_layer == 0){

#if 1   // local vars

            int     i_face_int, i_face_ext, i_dim,i_face, i_thread, i_cell, i_fluent_face;

            int     c0, c1;

            double  A[3];

            int     int_face_count = rcfd_get_int_face_count(handle, i_layer);

            int     *face_c0    = rcfd_get_face_c0  (handle, i_layer); // [i_face]

            int     *face_c1    = rcfd_get_face_c1  (handle, i_layer); // [i_face]

            double  **face_area = rcfd_get_face_area(handle, i_layer); // [i_face][i_dim]

            Domain  *d = Get_Domain(1);

            Thread  *t = NULL, *ct0 = NULL, *ct1 = NULL;
#endif

            // set _F(i_layer).c0/c1 and _F(i_layer).A
            {
                i_face_int = 0;

                i_face_ext = 0;

                thread_loop_f(t,d){if(THREAD_TYPE(t)==THREAD_F_INTERIOR){begin_f_loop(i_fluent_face, t){

                    c0 = (int)F_C0(i_fluent_face, t);

                    c1 = (int)F_C1(i_fluent_face, t);

                    ct0 = F_C0_THREAD(i_fluent_face, t);

                    ct1 = F_C1_THREAD(i_fluent_face, t);

                    if((c0 >= 0) && (c1 >= 0)){

                        if((C_PART(c0, ct0) == myid) && (C_PART(c1, ct1) == myid)){

                            i_face = i_face_int;

                            i_thread = Fluent_Topo->global_cell_thread_id[THREAD_ID(ct0)];

                            i_cell = Fluent_Topo->global_cell_id[i_thread][c0];

                            face_c0[i_face] = i_cell;

                            i_thread = Fluent_Topo->global_cell_thread_id[THREAD_ID(ct1)];

                            i_cell = Fluent_Topo->global_cell_id[i_thread][c1];

                            face_c1[i_face] = i_cell;

                            F_AREA(A, i_fluent_face, t);

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

                                face_area[i_face][i_dim] = A[i_dim];
                            }

                            i_face_int++;

                        }
                        else if((C_PART(c0, ct0) == myid) || (C_PART(c1, ct1) == myid)){

                            i_face = int_face_count + i_face_ext;

                            i_thread = Fluent_Topo->global_cell_thread_id[THREAD_ID(ct0)];

                            i_cell = Fluent_Topo->global_cell_id[i_thread][c0];

                            face_c0[i_face] = i_cell;

                            i_thread = Fluent_Topo->global_cell_thread_id[THREAD_ID(ct1)];

                            i_cell = Fluent_Topo->global_cell_id[i_thread][c1];

                            face_c1[i_face] = i_cell;

                            // eventually swap c0/c1 such that c1 is always corona cell
                            if(C_PART(c0, ct0) != myid){

                                i_cell = face_c0[i_face];

                                face_c0[i_face] = face_c1[i_face];

                                face_c1[i_face] = i_cell;
                            }

                            F_AREA(A, i_fluent_face, t);

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

                                face_area[i_face][i_dim] = A[i_dim];
                            }

                            i_face_ext++;
                        }
                    }

                }end_f_loop(i_fluent_face, t)}}
            }
        }
#endif // RP_NODE
    }

    /**
     * @brief extract the cell host information from the ANSYS Fluent mesh to
     *        initialize the corresponding rCFD cell data structures. Note that
     *        this overwrites the values that may have been read from rCFD topo
     *        files
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_init_cell_hosts(void *handle, const short i_layer)
    {
        if(i_layer == 0){    // only implemented for base layer

#if 1   // local vars

            int     i_cell, i_fluent_cell, i_thread;

            int     *cell_host = rcfd_get_cell_host(handle, i_layer);

            Domain  *d = Get_Domain(1);

            Thread  *t;
#endif

            thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){

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

                begin_c_loop_all(i_fluent_cell,t){

                    i_cell = Fluent_Topo->global_cell_id[i_thread][i_fluent_cell];

                    cell_host[i_cell] = C_PART(i_fluent_cell,t);

                }end_c_loop_all(i_fluent_cell,t)

            }}
        }
    }

    /**
     * @brief Convert ANSYS Fluent cell type and order of node indices to VTK definition
     * @param fluent_cell_type  Ansys Fluent cell type
     * @param cell_node_ids     Ansys Fluent cell node ids to be reordered to VTK definition
     * @return corresponding VTK cell type or vertex type if unsupported
     */
    int fluent_to_vtk_cell_type(int fluent_cell_type, int *cell_node_ids)
    {
        int tmp_node_id;
        switch(fluent_cell_type){

        case TRI_CELL:
            // Fluent = VTK
            return RCFD_CELL_TRIANGLE;
#if RP_2D
        case QUAD_CELL:
            // Fluent = VTK
            return RCFD_CELL_QUADRILATERAL;
#endif
        case TET_CELL:
            // Fluent =/= VTK
            tmp_node_id = cell_node_ids[2];
            cell_node_ids[2] = cell_node_ids[3];
            cell_node_ids[3] = tmp_node_id;
            return RCFD_CELL_TETRAHEDRON;
        case PYRAMID_CELL:
            // Fluent =/= VTK
            tmp_node_id = cell_node_ids[2];
            cell_node_ids[2] = cell_node_ids[4];
            cell_node_ids[4] = tmp_node_id;
            return RCFD_CELL_PYRAMID;
        case WEDGE_CELL:
            // Fluent = VTK
            return RCFD_CELL_WEDGE;
        case HEX_CELL:
            tmp_node_id = cell_node_ids[2];
            cell_node_ids[2] = cell_node_ids[4];
            cell_node_ids[4] = tmp_node_id;
            tmp_node_id = cell_node_ids[6];
            cell_node_ids[6] = cell_node_ids[8];
            cell_node_ids[8] = tmp_node_id;
            return RCFD_CELL_HEXAHEDRON;
        case POLY_CELL:
            return RCFD_CELL_POLYHEDRON;
        default:
            // store any unsupported cell type as vertex
            Message("\n... WARNING: fluent_to_vtk_cell_type: unsupported cell type %d", fluent_cell_type);
            return RCFD_CELL_VERTEX;
        }
    }

    /**
     * @brief extract the cell type and cell nodes from the ANSYS Fluent mesh
     *        to initialize the corresponding rCFD cell/node data structures.
     *
     * @param handle    pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_init_nodes(void *handle, const short i_layer)
    {
#if RP_NODE
        if(i_layer == 0){

            int i_node, i_fluent_cell, i_thread;
            int node_local_id = 0;
            double  **node_x  = rcfd_get_node(handle, i_layer);
            Node    *node;
            Domain  *d = Get_Domain(1);
            Thread  *t = NULL;

            if(!node_x){

                Message("\n... WARNING: [%d] rCFD has not allocated memory for cell nodes\n", myid);
                return;
            }

            // init node properties for all nodes
            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

                begin_c_loop_all(i_fluent_cell, t){

                    c_node_loop(i_fluent_cell,t,i_node){

                        node = C_NODE(i_fluent_cell,t,i_node);

                        N_UDMI(node, NODE_VSITED)   =  0.0; // not yet visited
                        N_UDMI(node, NODE_LOCAL_ID) = -1.0; // not yet set
                    }
                }end_c_loop_all(i_fluent_cell,t)
            }}

            // set node positions and local node IDs
            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

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

                begin_c_loop_all(i_fluent_cell, t){

                    c_node_loop(i_fluent_cell,t,i_node){

                        node = C_NODE(i_fluent_cell,t,i_node);

                        if(N_UDMI(node, NODE_VSITED) == 0.0){

                            node_x[node_local_id][0] = NODE_X(node);
                            node_x[node_local_id][1] = NODE_Y(node);
                            node_x[node_local_id][2] = NODE_Z(node);

                            N_UDMI(node, NODE_VSITED)   = 1.0;
                            N_UDMI(node, NODE_LOCAL_ID) = node_local_id;
                            node_local_id++;
                        }
                    }
                }end_c_loop_all(i_fluent_cell,t)
            }}

            // set cell types and node indices
            thread_loop_c(t, d){if(FLUID_CELL_THREAD_P(t)){

                int i_cell, i_face, cell_type;
                int number_of_cell_faces;

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

                begin_c_loop_all(i_fluent_cell, t){

                    int *poly_face_nodes = NULL; // contains number of poly faces followed by number of face nodes
                    int *cell_node_ids = NULL; // contains number of nodes followed by node IDs

                    i_cell = Fluent_Topo->global_cell_id[i_thread][i_fluent_cell];

                    cell_type = C_TYPE(i_fluent_cell, t);

                    number_of_cell_faces = C_NFACES(i_fluent_cell, t); // cell face count

                    if(number_of_cell_faces > RCFD_MAX_FACES_PER_POLY_CELL){

                        Message("\n... WARNING: [%d] rCFD supports max. %d faces per polyhedron, %d were detected in mesh\n",
                                myid, RCFD_MAX_FACES_PER_POLY_CELL, number_of_cell_faces);
                    }

                    if(cell_type == POLY_CELL){

                        poly_face_nodes = malloc(sizeof(int) * (1 + number_of_cell_faces));
                        poly_face_nodes[0] = number_of_cell_faces;
                        // Note that C_NNODES(i_fluent_cell, t) is different from the sum of all face nodes F_NNODES(i_fluent_face, tf)
                        cell_node_ids = (int*)malloc(sizeof(int) * 1);
                        cell_node_ids[0] = 0;

                        c_face_loop(i_fluent_cell,t,i_face){

                            int i_fluent_face = C_FACE(i_fluent_cell, t, i_face);
                            Thread *tf = C_FACE_THREAD(i_fluent_cell, t, i_face);
                            int number_of_face_nodes = F_NNODES(i_fluent_face, tf);

                            poly_face_nodes[i_face+1] = number_of_face_nodes;
                            cell_node_ids = realloc(cell_node_ids, sizeof(int) * (1 + cell_node_ids[0] + number_of_face_nodes));

                            f_node_loop(i_fluent_face, tf, i_node){

                                node = F_NODE(i_fluent_face, tf, i_node);
                                cell_node_ids[1+cell_node_ids[0]+i_node] = N_UDMI(node, NODE_LOCAL_ID);
                            }
                            cell_node_ids[0] += number_of_face_nodes;
                        }
                    }
                    else{

                        cell_node_ids = malloc(sizeof(int) * (1 + C_NNODES(i_fluent_cell, t)));
                        cell_node_ids[0] = C_NNODES(i_fluent_cell, t); // cell node count

                        c_node_loop(i_fluent_cell,t,i_node){

                            if(i_node == RCFD_MAX_NODES_PER_CELL) break;

                            node = C_NODE(i_fluent_cell,t,i_node);

                            cell_node_ids[i_node+1] = N_UDMI(node, NODE_LOCAL_ID);
                        }
                    }

                    cell_type = fluent_to_vtk_cell_type(cell_type, cell_node_ids);

                    rcfd_set_cell_shape(handle, i_layer, i_cell, cell_type, cell_node_ids, poly_face_nodes);

                    free(poly_face_nodes);
                    free(cell_node_ids);

                }end_c_loop_all(i_fluent_cell,t)
            }}
        }

#endif // RP_NODE
    }

    /**
     * @brief get the ANSYS Fluent phase ID from the given rCFD phase index
     *
     * @param rcfd      pointer to previously created rCFD instance
     * @param i_phase   specifies the rCFD phase index
     * @return the ANSYS Fluent phase ID
     */
    short fluent_get_phase_id(rCFD_type *rcfd, const short i_phase)
    {
#if RP_NODE
        if(rcfd->Solver_Dict->single_phase_vof){

            return rcfd->Solver_Dict->single_phase_vof_id;
        }
        else{

            return i_phase;
        }
#endif
        return 0;
    }

    /**
     * @brief extract the cell VOF information from the ANSYS Fluent mesh to
     *        initialize the corresponding rCFD VOF data structures in
     *        preparation mode
     *
     * @param rcfd      pointer to previously created rCFD instance
     * @param i_layer   specifies the cell layer
     */
    void fluent_set_vof(rCFD_type *rcfd, const short i_layer)
    {
#if RP_NODE

#if 1
        int     i_cell, i_fluent_cell, i_thread, i_phase, i_phase_fluent;
        int     phase_count = rcfd_get_phase_count(myrcfd);
        double  **cell_vof = rcfd_get_cell_vof(myrcfd, i_layer, 0);

        Domain  *d = Get_Domain(1);
        Thread  *t = NULL, *t_phase = NULL;
#endif

        if(rcfd->Solver_Dict->mode == RCFD_MODE_RUN){

            rcfd_report0(myrcfd, RCFD_REPORT_WARNING,
                         "fluent_set_vof: only sets first frame's vof (not designed for run mode)");
            return;
        }

        thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){

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

            begin_c_loop_all(i_fluent_cell,t){

                i_cell = Fluent_Topo->global_cell_id[i_thread][i_fluent_cell];

                for(i_phase = 0; i_phase < phase_count; i_phase++){

                    i_phase_fluent = fluent_get_phase_id(rcfd, i_phase);

                    t_phase = THREAD_SUB_THREAD(t, i_phase_fluent);

                    cell_vof[i_cell][i_phase] = C_VOF(i_fluent_cell, t_phase);
                }
            }end_c_loop_all(i_fluent_cell,t)
        }}
#endif
    }

#endif // RCFD_FLUENT_H
