#define BUF_LIM 201
#define VAR_LIM 31
#define NODES_PER_ELEMENT 4

#define LOAD_ALL 0
#define LOAD_REDUCED 1
#define LOAD_GLOBAL 2

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <conio.h>
#include <vector>
#include <algorithm>

struct variables
    {
        FILE * F = NULL;   // variable declarations
        char * fil_buf = NULL;
        char * name_buf = NULL;
        char * var_buf = NULL;

        FILE * S = NULL;    // substructure bin info
        FILE * R = NULL;    // register

        double E = 0;
        double nu = 0;
        double force_multiplier = 0;    // constant for making forces smaller on paper but just right in calculations
        double density = 0;
        double thickness = 0;
        double scale_factor = 0;
        double ** C = NULL;    // the modulus matrix from 'sigma = C * epsilon'

        unsigned long number_of_nodes = 0;  // total number of nodes of the simulation
        unsigned long number_of_elements = 0;   // total number of elements of the simulation
        unsigned long number_of_substructures = 0;  // total number of substructures of the simulation

        unsigned long number_of_support_types = 0;  // total number of different types of supports
        unsigned long ** support_nodes = NULL; // indices of nodes that belong to a given support type
        unsigned long * support_index = NULL;  // array that stores indices of support type in case the model and program indication differ
        char ** support_names = NULL;  // array for storing the names of given support types
        unsigned long * number_of_support_nodes_per_type = NULL;    // array pointer that stores the lengths of different support type arrays of nodes
        unsigned long ** element_nodes = NULL; // stores the indices of nodes that belong to given elements
        unsigned long * number_of_substructure_elements = NULL;    // stores the amount of elements in given substructures
        unsigned long * substructure_elements = NULL;     // stores indices of elements that belong to given substructures

        unsigned long number_of_forces = 0; // the number of forces involved in the simulation
        unsigned long * force_indices = NULL;    // stores the indices of nodes where forces act
        double ** forces = NULL;    // stores the vectors of forces acting on nodes

        unsigned long global_dof = 0;

        double ** stiffness = NULL;
        double ** stiffness_reduced = NULL;
        double * displacements = NULL;
        double ** mass = NULL;
        double * node_coordinates_x = NULL;
        double * node_coordinates_y = NULL;
        double * force = NULL;
        double * force_reduced = NULL;
        double ** reactions = NULL;

        double ** locations = NULL;
        double ** weights = NULL;

        double *** stresses = NULL;

        unsigned long fixed_index = 0;

        unsigned long number_H = 0;
        unsigned long number_V = 0;
        unsigned long * H = NULL;
        unsigned long * V = NULL;
        unsigned long * sub_H = NULL;
        unsigned long * sub_V = NULL;
        unsigned long bound_length = 0;
        unsigned long * substructure_nodes = NULL;
        unsigned long ** substructure_element_nodes = NULL;
        unsigned long * substructure_supports = NULL;
        unsigned long number_of_substructure_support_nodes = 0;

        double *** g_stiffness = NULL;
        double *** g_stiffness_reduced = NULL;
        double ** g_force = NULL;
        double ** g_force_reduced = NULL;
        unsigned long * g_num_H = NULL;
        unsigned long * g_num_V = NULL;
        unsigned long ** g_H = NULL;
        unsigned long ** g_V = NULL;

        unsigned long * red_nodes_H = NULL;
        unsigned long num_of_red_nodes_H = 0;
        double ** g_K_h = NULL;
        double * g_F_h = NULL;
        unsigned long * g_prescribed_dof = NULL;
        double ** g_U_h = NULL;
        double * g_U_v = NULL;

/////////////// utility variables

        unsigned long num_buf = 0;  // used to buffer numerical values for comparison purposes
        int breaker = 0;    // used in sorting to indicate a match
        unsigned long counter = 0;

        const char * file_substructures = "matrix_storage\\substructures.bin";
        const char * file_register = "matrix_storage\\register.bin";
        const char * file_inverse = "matrix_storage\\inverse\\inverse_";
        const char * file_bin_ext = ".bin";
    };

#include "Head.h"

int main()
{
    variables I;
    time_t ti;
    time(&ti);
    FILE * L = NULL;
    char * logname = NULL;
    logname = (char *)calloc(strlen("logs/log.txt"), sizeof(char));

    strcpy(logname, "logs/log.txt");

    L = fopen(logname, "w");

    if(L == NULL)
    {
        printf("Failed to open log file.\n%s.\n", logname);
    }
    free(logname);

    fprintf(L, "%s************************\n\n", ctime(&ti));

// // // // // // // // // // // // // // // // // // // // // // // //// // // // // // // // // // // // // // // // // // // // // // // //

    I.F = fopen("setup.txt", "r");
    if(I.F == NULL)   // if the file does not exist, terminte the program and generate new setup file
    {
        printf("setup.txt not found.\nGenerating the file for you.\nYou need to configure the values before\nproceeding to starting this simulation again.\n");
        setup_gen();    // create new setup file
        return -1;  // terminate the program
    }

    I.fil_buf = (char*)calloc(BUF_LIM, sizeof(char)); // allocate BUF_LIM amount of memory for reading the input files
    while(!feof(I.F)) // run until the end of the file
    {
        fgets(I.fil_buf, BUF_LIM - 1, I.F); // filling the reading buffer with a line from the file
        I.name_buf = parse_input(I.fil_buf, '<', '>');  // finding the variable
        I.var_buf = parse_input(I.fil_buf, '{', '}');   // and its value

        if(strncmp(I.name_buf, "E", 1) == 0)    // if E is found
        {
            I.E = strtod(I.var_buf, NULL);  // converts the written value to double precision floating point
        }else if(strncmp(I.name_buf, "nu", 2) == 0)    // if it found nu
        {
            I.nu = strtod(I.var_buf, NULL); // converts the written value to double precision floating point
        }else if(strncmp(I.name_buf, "force_multiplier", 16) == 0)
        {
            I.force_multiplier = strtod(I.var_buf, NULL);
        }else if(strncmp(I.name_buf, "density", 7) == 0)
        {
            I.density = strtod(I.var_buf, NULL);
        }else if(strncmp(I.name_buf, "thickness", 9) == 0)
        {
            I.thickness = strtod(I.var_buf, NULL);
        }else if(strncmp(I.name_buf, "scale_factor", 12) == 0)
        {
            I.scale_factor = strtod(I.var_buf, NULL);
        }

        memset(I.fil_buf, 0, BUF_LIM);    // clear the file reading buffer
        free(I.name_buf); // deallocates buffers so that they can be allocated in the next iteration
        free(I.var_buf);
    }
    fclose(I.F);  // closes the setup.txt

    if(I.E == 0 || I.nu == 0 || I.force_multiplier == 0 || I.density == 0 || I.thickness == 0 || I.scale_factor == 0)   // final check of the values read from setup.txt
    {
        printf("Failed to read the setup.txt file.\nE = %f;\nnu = %f;\nforce_multiplier = %f;\ndensity = %f;\nthickness = %f;\nscale_factor = %f;\n", I.E, I.nu, I.force_multiplier, I.density, I.thickness, I.scale_factor);   // error message
        return -1;
    }

    I.C = makeC(I.E, I.nu);

// // // // // // // // // // // // // // // // // // // // // // // //// // // // // // // // // // // // // // // // // // // // // // // //

    I.F = fopen("model_superelementary.mkp", "r");    // opening the model.mkp file with geometry info
    if(I.F == NULL)   // if it did not open
    {
        printf("model_superelementary.mkp not found.\nMake sure the file is located in the same folder as this program's executable.\n");    // error message
        return -2;
    }

    while(!feof(I.F)) // if the end of the file was not reached yet
    {
        fgets(I.fil_buf, BUF_LIM - 1, I.F); // read a line from the file and store it in the buffer

        if(I.fil_buf[0] == '%' || strlen(I.fil_buf) == 1)   // if the read line is a comment or only contains newline character (aka is empty)
        {
            continue;   // read another line
        }else   // if the line contains relevant info...
        {
            if(strncmp(I.fil_buf, ".<number_of_nodes>", 18) == 0)  // if the line contains the header for the node info
            {
                I.var_buf = parse_input(I.fil_buf, '{', '}');   // get the number of nodes
                I.number_of_nodes = strtoul(I.var_buf, NULL, 10);    // extract the number of nodes and store it in corresponding variable

                free(I.var_buf);  // free the line reading buffer for further use


            }else if(strncmp(I.fil_buf, "#<number_of_elements>", 21) == 0) // if the line contains the header for the element info
            {
                I.var_buf = parse_input(I.fil_buf, '{', '}');   // get the number of elements
                I.number_of_elements = strtoul(I.var_buf, NULL, 10); // extract the number of elements and store it in corresponding variable

                free(I.var_buf);  // free the line reading buffer for further use

            }else if(strncmp(I.fil_buf, "_<number_of_substructures>", 26) == 0)
            {
                I.var_buf = parse_input(I.fil_buf, '{', '}');   // if the line contains substructure info
                I.number_of_substructures = strtoul(I.var_buf, NULL, 10);    // extract the number of substructures and store it in corresponding variable

                free(I.var_buf);

            }else if(strncmp(I.fil_buf, "^<number_of_support_types>", 26) == 0)
            {
                I.var_buf = parse_input(I.fil_buf, '{', '}');   // if the line contains support info
                I.number_of_support_types = strtoul(I.var_buf, NULL, 10);    // extract the number of support types and store it in corresponding variable

                free(I.var_buf);
                memset(I.fil_buf, 0, BUF_LIM);    //clear the file reading buffer

                I.support_nodes = (unsigned long **)calloc(I.number_of_support_types, sizeof(unsigned long));  // allocate space for support nodes indices
                I.support_index = (unsigned long *)calloc(I.number_of_support_types, sizeof(unsigned long)); // in case model file numbers support types in different order or a whole new way
                I.support_names = (char **)calloc(I.number_of_support_types, sizeof(char)); // allocate array for storing names of the support types
                I.number_of_support_nodes_per_type = (unsigned long *)calloc(I.number_of_support_types, sizeof(unsigned long));    // allocate array for storing amounts of nodes in different support types

                for(unsigned long i = 0; i < I.number_of_support_types; i++)  // loop for collecting node indices of given support types
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F); // read a line from the file and store it in the buffer

                    if(strncmp(I.fil_buf, "^\"", 1) == 0) // if the line contains support info
                    {
                        I.var_buf = parse_input(I.fil_buf, '"', '\'');  // find the name of the support type...
                        I.support_names[i] = (char*)calloc(strlen(I.var_buf), sizeof(char));    // ...get its length...
                        strcpy(I.support_names[i], I.var_buf);  // ...and store it

                        free(I.var_buf);  // clear the buffer for further use

                        I.var_buf = parse_input(I.fil_buf, '<', '>');   // find the index of the support type...
                        I.support_index[i] = strtoul(I.var_buf, NULL, 10);  // ...and store it

                        free(I.var_buf);  // clear the buffer for further use

                        I.var_buf = parse_input(I.fil_buf, '{', '}');   // get the amount of nodes that at as supports of the given type...
                        I.number_of_support_nodes_per_type[i] = strtoul(I.var_buf, NULL, 10);   // ...and store it

                        free(I.var_buf);  // free the buffer for further use
                        memset(I.fil_buf, 0, BUF_LIM);    // clear the file reading buffer as we know everything from current line

                        I.support_nodes[i] = (unsigned long *)calloc(I.number_of_support_nodes_per_type[i], sizeof(unsigned long));    // allocate an appropriate amount of memory for the support type node array

                        for(unsigned long j = 0; j < I.number_of_support_nodes_per_type[i]; j++)  // looping through nodes in given support type
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F); // read another line

                            I.var_buf = parse_input(I.fil_buf, '^', '=');   // look for the index of support type...
                            I.num_buf = strtoul(I.var_buf, NULL, 10);    // ...and store it temporarily

                            free(I.var_buf);  // clear the buffer for further use

                            I.var_buf = parse_input(I.fil_buf, '{', '}');   // look for the the node index
                            I.support_nodes[i][j] = strtoul(I.var_buf, NULL, 10);   // store the index number in corresponding array

                            free(I.var_buf);

                            if(I.num_buf != I.support_index[i]) // check if the index makes sense
                            {
                                printf("Support index %lu in its %luth node.\n", I.support_index[i], I.support_nodes[i][j]); // if the index does not match, print error and terminate
                                return -3;
                            }

                            memset(I.fil_buf, 0, BUF_LIM);    // clear the file reading buffer so that other processes get an empty buffer
                        }
                    }
                    memset(I.fil_buf, 0, BUF_LIM);    // clear the file reading buffer
                }
            }else if(strncmp(I.fil_buf, "$<number_of_forces>", 19) == 0)
            {
                I.var_buf = parse_input(I.fil_buf, '{', '}');
                I.number_of_forces = strtoul(I.var_buf, NULL, 10);

                free(I.var_buf);

                I.force_indices = (unsigned long *)calloc(I.number_of_forces, sizeof(unsigned long));
                I.forces = (double **)calloc(I.number_of_forces, sizeof(double));

                for(unsigned long i = 0; i < I.number_of_forces; i++)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);

                    I.force_indices[i] = strtoul(parse_input(I.fil_buf, '$', '='), NULL, 10);

                    I.forces[i] = (double *)calloc(2, sizeof(double));
                    I.forces[i][0] = strtod(parse_input(I.fil_buf, '{', ';'), NULL);   // storing the x force component
                    I.forces[i][1] = strtod(parse_input(I.fil_buf, ';', '}'), NULL);   // storing the y force component
                }
            }else if(strncmp(I.fil_buf, "!*", 2) == 0)
            {
                printf("Input complete.\n");
            }
        }

        memset(I.fil_buf, 0, BUF_LIM);    //clear the file reading buffer
    }

    for(unsigned long i = 0; i < I.number_of_support_types; i++)
    {
        if((strcmp(I.support_names[i], "fixed") == 0) || (strcmp(I.support_names[i], "Fixed") == 0))    // finding the index of fixed support
        {
            I.fixed_index = i;
            printf("fixed found %lu\n", I.fixed_index);
            break;
        }
    }

    I.locations = gaussQuadratureLoc();
    I.weights = gaussQuadratureWei();
    I.number_of_substructure_elements = (unsigned long *)calloc(I.number_of_substructures, sizeof(unsigned long));

/// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    printf("Proceeding to the superelementary part\n");

    I.S = fopen(I.file_substructures, "w+b");
    if(I.S == NULL)
    {
        printf("Failed to open %s\nMaybe the folder is missing.", I.file_substructures);
        return -1;
    }
    I.R = fopen(I.file_register, "w+b");
    if(I.R == NULL)
    {
        printf("Failed to open %s\nMaybe the folder is missing.", I.file_register);
        fclose(I.S);
        return -1;
    }

    for(unsigned long cur_sub = 0; cur_sub < I.number_of_substructures; cur_sub++)    // all the substructures
    {
        fseek(I.F, 0, SEEK_SET);
        while(true)
        {
            fgets(I.fil_buf, BUF_LIM - 1, I.F); // read a line from the file and store it in the buffer

            if(strncmp(I.fil_buf, "_<number_of_substructures>", 26) == 0)  // if the line contains the header for the substructure element info
            {
                while(true)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    if(strncmp(I.fil_buf, "_<", 2) == 0)
                    {
                        I.var_buf = parse_input(I.fil_buf, '<', '>');
                        if(strtoul(I.var_buf, NULL, 10) == cur_sub)   // if the file is beginning the current substructure block
                        {
                            printf("found substructure %lu\n", cur_sub);
                            I.var_buf = parse_input(I.fil_buf, '{', '}');
                            I.number_of_substructure_elements[cur_sub] = strtoul(I.var_buf, NULL, 10);
                            I.substructure_elements = (unsigned long *)calloc(I.number_of_substructure_elements[cur_sub], sizeof(unsigned long));

                            for(unsigned long e = 0; e < I.number_of_substructure_elements[cur_sub]; e++)
                            {
                                fgets(I.fil_buf, BUF_LIM - 1, I.F);
                                I.substructure_elements[e] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);
                            }

                            if(I.fil_buf[strlen(I.fil_buf) - 2] == ';')
                            {
                                printf("elements from sub %lu OK\n", cur_sub);
                            }else
                            {
                                printf("elements from sub %lu not quite there...\n", cur_sub);
                                printf("%s\n", I.fil_buf);
                            }
                            break;
                        }
                    }
                }

                fseek(I.F, 0, SEEK_SET);
                while(true)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    if(strncmp(I.fil_buf, "#<number_of_elements>", 21) == 0)  // if the line contains the header for the element info
                    {
                        I.element_nodes = (unsigned long **)calloc(I.number_of_substructure_elements[cur_sub], sizeof(unsigned long));
                        I.counter = 0;
                        while(I.counter < I.number_of_substructure_elements[cur_sub])
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F);
                            I.var_buf = parse_input(I.fil_buf, '#', '=');
                            if(strtoul(I.var_buf, NULL, 10) == I.substructure_elements[I.counter])
                            {
                                I.element_nodes[I.counter] = (unsigned long *)calloc(NODES_PER_ELEMENT, sizeof(unsigned long));
                                I.element_nodes[I.counter][0] = strtoul(parse_input(I.fil_buf, '{', ','), NULL, 10);
                                I.element_nodes[I.counter][1] = strtoul(parse_input(I.fil_buf, ',', ':'), NULL, 10);
                                I.element_nodes[I.counter][2] = strtoul(parse_input(I.fil_buf, ':', ';'), NULL, 10);
                                I.element_nodes[I.counter][3] = strtoul(parse_input(I.fil_buf, ';', '}'), NULL, 10);
                                I.counter++;
                            }
                        }
                        break;
                    }
                }

                while(true)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    if(strncmp(I.fil_buf, "&<substructure_nodes>", 21) == 0)
                    {
                        while(true)
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F);
                            if(strtoul(parse_input(I.fil_buf, '<', '>'), NULL, 10) == cur_sub)
                            {
                                I.num_buf = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);
                                I.substructure_nodes = (unsigned long *)calloc(I.num_buf, sizeof(unsigned long));

                                for(unsigned long x = 0; x < I.num_buf; x++)
                                {
                                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                                    I.substructure_nodes[x] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);
                                }
                                break;
                            }
                        }
                        break;
                    }
                }

                fseek(I.F, 0, SEEK_SET);
                while(true)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    if(strncmp(I.fil_buf, "|<substructure_nodes_bounded>", 29) == 0)
                    {
                        while(true)
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F);
                            if(strtoul(parse_input(I.fil_buf, '<', '>'), NULL, 10) == cur_sub)
                            {
                                break;
                            }
                        }

                        I.number_V = strtoul(parse_input(I.fil_buf, '[', ']'), NULL, 10);
                        I.number_H = strtoul(parse_input(I.fil_buf, '(', ')'), NULL, 10);
                        I.V = (unsigned long *)calloc(I.number_V, sizeof(unsigned long));
                        for(unsigned long v = 0; v < I.number_V; v++)
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F);
                            I.V[v] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);
                        }

                        I.H = (unsigned long *)calloc(I.number_H, sizeof(unsigned long));
                        for(unsigned long h = 0; h < I.number_H; h++)
                        {
                            fgets(I.fil_buf, BUF_LIM - 1, I.F);
                            I.H[h] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);
                        }
                        break;
                    }
                }

                fseek(I.F, 0, SEEK_SET);
                while(true)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    if(strncmp(I.fil_buf, ".<number_of_nodes>", 18) == 0)
                    {
                        I.node_coordinates_x = (double *)calloc(I.number_V + I.number_H, sizeof(double));
                        I.node_coordinates_y = (double *)calloc(I.number_V + I.number_H, sizeof(double));
                        I.counter = 0;
                        break;
                    }
                }
                while(I.counter < I.number_V + I.number_H)
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);
                    I.num_buf = strtoul(parse_input(I.fil_buf, '.', '='), NULL, 10);

                    if(is_in(I.num_buf, I.V, I.number_V) || is_in(I.num_buf, I.H, I.number_H))
                    {
                        I.node_coordinates_x[I.counter] = strtod(parse_input(I.fil_buf, '{', ';'), NULL);
                        I.node_coordinates_y[I.counter] = strtod(parse_input(I.fil_buf, ';', '}'), NULL);
                        I.counter++;
                    }
                }
                break;
            }
        }

        I.number_of_substructure_support_nodes = sub_supports(I.substructure_supports, I.support_nodes[I.fixed_index], I.substructure_nodes, I.number_of_support_nodes_per_type[I.fixed_index], I.number_H + I.number_V);
        I.substructure_element_nodes = sub_nodes(I.element_nodes, I.substructure_nodes, I.number_of_substructure_elements[cur_sub], I.number_H + I.number_V);
        sub_bounds(I.sub_H, I.sub_V, I.substructure_nodes, I.H, I.V, I.number_H, I.number_V);

/// the end of planned per structure section
/// begin the reduced stiffness calculation

        I.global_dof = 2 * (I.number_H + I.number_V);
        formStiffness2D(I.stiffness, I.mass, I.global_dof, I.number_of_substructure_elements[cur_sub], I.substructure_element_nodes, I.number_H + I.number_V, I.node_coordinates_x, I.node_coordinates_y, I.C, I.density, I.thickness, I.locations, I.weights, L);
        I.force = sub_force(I.force_indices, I.forces, I.substructure_nodes, I.global_dof, I.number_H + I.number_V, I.number_of_forces, I.force_multiplier);
        reorder(I.stiffness, I.force, I.substructure_element_nodes[cur_sub], I.sub_V, I.sub_H, I.number_V, I.number_H);
        I.stiffness_reduced = matCopy(I.stiffness, 2 * (I.number_H + I.number_V), 2 * (I.number_H + I.number_V));
        I.force_reduced = vecCopy(I.force, 2 * (I.number_H + I.number_V));
        reduce_H(I.stiffness_reduced, I.force_reduced, I.number_H, I.number_V, cur_sub, I.file_inverse, I.file_bin_ext, L);

/// store the reduced matrices and all necessary substructure info

        store_reduced(I.S, I.R, I.stiffness, I.force, I.stiffness_reduced, I.force_reduced, I.number_H, I.number_V, I.H, I.V);

/// get rid of all the other stuff

        for(unsigned long i = 0; i < 2 * I.number_H; i++)
        {
            free(I.stiffness_reduced);
        }

        for(unsigned long i = 0; i < I.number_of_substructure_elements[cur_sub]; i++)
        {
            free(I.element_nodes[i]);
        }

        for(unsigned long i = 0; i < I.global_dof; i++)
        {
            free(I.stiffness[i]);
            free(I.mass[i]);
        }

        free(I.stiffness_reduced);
        free(I.force_reduced);
        free(I.stiffness);
        free(I.mass);
        free(I.node_coordinates_x);
        free(I.node_coordinates_y);
        free(I.H);
        free(I.V);
        free(I.sub_H);
        free(I.sub_V);
        free(I.element_nodes);
        free(I.substructure_elements);
        free(I.substructure_nodes);
        free(I.substructure_element_nodes);
        free(I.force);
        free(I.substructure_supports);
    }
    fclose(I.F);
    fclose(I.S);
    fclose(I.R);

    for(unsigned long i = 0; i < I.number_of_forces; i++)
    {
        free(I.forces[i]);
    }
    for(unsigned long i = 0; i < I.number_of_support_types; i++)
    {
        free(I.support_names);
    }
    free(I.support_names);
    free(I.support_index);
    free(I.forces);
    free(I.force_indices);

/// calculate the global reduced stiffness matrix
    printf("Proceeding to calculating the reduced part\n");

    I.S = fopen(I.file_substructures, "rb");
    if(I.S == NULL)
    {
        printf("Failed to open %s for READING...\n", I.file_substructures);
        return -1;
    }
    I.R = fopen(I.file_register, "rb");
    if(I.R == NULL)
    {
        printf("Failed to open %s for READING...\n", I.file_register);
        return -1;
    }

    I.g_stiffness_reduced = (double ***)calloc(I.number_of_substructures, sizeof(double));
    I.g_force_reduced = (double **)calloc(I.number_of_substructures, sizeof(double));
    I.g_H = (unsigned long **)calloc(I.number_of_substructures, sizeof(unsigned long));
    I.g_V = (unsigned long **)calloc(I.number_of_substructures, sizeof(unsigned long));
    I.g_num_H = (unsigned long *)calloc(I.number_of_substructures, sizeof(unsigned long));
    I.g_num_V = (unsigned long *)calloc(I.number_of_substructures, sizeof(unsigned long));

    for(unsigned long i = 0; i < I.number_of_substructures; i++)
    {
        load_reduced(I.S, I.R, i, I.stiffness, I.force, I.g_stiffness_reduced[i], I.g_force_reduced[i], I.g_num_H[i], I.g_num_V[i], I.g_H[i], I.g_V[i], LOAD_REDUCED);
    }

/// generate global K reduced and calculate U_h

    globalize_H(I.red_nodes_H, I.num_of_red_nodes_H, I.g_H, I.g_num_H, I.number_of_substructures);
    global_K_h(I.g_K_h, I.red_nodes_H, I.num_of_red_nodes_H, I.g_stiffness_reduced, I.g_H, I.g_num_H, I.number_of_substructures);
    global_prescribed_dof(I.g_prescribed_dof, I.support_nodes[I.fixed_index], I.number_of_support_nodes_per_type[I.fixed_index], I.red_nodes_H, I.num_of_red_nodes_H);
    global_F_h(I.g_F_h, I.red_nodes_H, I.num_of_red_nodes_H, I.g_force_reduced, I.g_H, I.g_num_H, I.number_of_substructures);

/// time for the grand calculation

    I.g_U_h = solution(2 * I.num_of_red_nodes_H, I.red_nodes_H, I.g_prescribed_dof, I.number_of_support_nodes_per_type[I.fixed_index], I.g_K_h, I.g_F_h, I.num_of_red_nodes_H,/*I.active_dof, I.num_active_dof, */L);
    I.displacements = (double *)calloc(2*I.number_of_nodes, sizeof(double));
    store_disp_H(I.g_U_h, I.displacements, I.red_nodes_H, I.num_of_red_nodes_H, I.number_of_nodes);

    for(unsigned long i = 0; i < 2*I.num_of_red_nodes_H; i++)
    {
        free(I.g_K_h[i]);
    }
    free(I.g_K_h);
    free(I.g_prescribed_dof);
    free(I.g_F_h);

    for(unsigned long i = 0; i < I.number_of_support_types; i++)
    {
        free(I.support_nodes);
    }
    free(I.support_nodes);

    for(unsigned long i = 0; i < I.number_of_substructures; i++)
    {

        for(unsigned long j = 0; j < 2 * (I.g_num_H[i]); j++)
        {
            free(I.g_stiffness_reduced[i][j]);
        }

        free(I.g_stiffness_reduced[i]);
        free(I.g_force_reduced[i]);
        free(I.g_H[i]);
        free(I.g_V[i]);
    }

    free(I.g_stiffness_reduced);
    free(I.g_force_reduced);
    free(I.g_H);
    free(I.g_V);

    free(I.g_num_H);
    free(I.g_num_V);

/// calculate individual U_vv's

    I.g_stiffness_reduced = (double ***)calloc(1, sizeof(double));
    I.g_force_reduced = (double **)calloc(1, sizeof(double));
    I.g_H = (unsigned long **)calloc(1, sizeof(unsigned long));
    I.g_V = (unsigned long **)calloc(1, sizeof(unsigned long));
    I.g_num_H = (unsigned long *)calloc(1, sizeof(unsigned long));
    I.g_num_V = (unsigned long *)calloc(1, sizeof(unsigned long));

    for(unsigned long s = 0; s < I.number_of_substructures; s++)
    {

        load_reduced(I.S, I.R, s, I.stiffness, I.force, I.g_stiffness_reduced[0], I.g_force_reduced[0], I.g_num_H[0], I.g_num_V[0], I.g_H[0], I.g_V[0], LOAD_GLOBAL);

        if(I.g_stiffness_reduced != NULL)
        {
            free(I.g_stiffness_reduced);
        }
        if(I.g_force_reduced != NULL)
        {
            free(I.g_force_reduced);
        }

/// calculate the individual U_vs

        reduce_V(I.g_U_v, I.stiffness, I.force, I.g_U_h, /*I.active_dof,*/ I.g_num_H[0], I.g_num_V[0], /*I.num_active_dof,*/ I.g_H[0], I.red_nodes_H, I.num_of_red_nodes_H, s, I.file_inverse, I.file_bin_ext, L);
        store_disp_V(I.g_U_v, I.displacements, I.g_V[0], I.g_num_V[0], I.number_of_nodes);

        for(unsigned long j = 0; j < 2 * (I.g_num_H[0] + I.g_num_V[0]); j++)
        {
            free(I.stiffness[j]);
        }
        free(I.stiffness);
        free(I.force);
        free(I.g_H[0]);
        free(I.g_V[0]);
        free(I.g_U_v);
    }

    I.F = fopen("model_superelementary.mkp", "r");    // opening the model.mkp file with geometry info
    if(I.F == NULL)   // if it did not open
    {
        printf("model_superelementary.mkp not found.\nMake sure the file is located in the same folder as this program's executable.\n");    // error message
        return -2;
    }

    rewind(I.F);
    load_stuff(I.node_coordinates_x, I.node_coordinates_y, I.element_nodes, I.F, I.number_of_nodes, I.number_of_elements);
    fclose(I.F);

    I.stresses = stresses2D(I.global_dof, I.number_of_elements, I.element_nodes, I.number_of_nodes, I.node_coordinates_x, I.node_coordinates_y, I.displacements, I.C);

    make_solutions(I.number_of_nodes, I.number_of_elements, I.element_nodes, I.node_coordinates_x, I.node_coordinates_y, I.displacements, I.stresses, I.scale_factor, ctime(&ti));

    for(unsigned long i = 0; i < 3; i++)
    {
        free(I.C[i]);
    }

    for(unsigned long i = 0; i < I.number_of_elements; i++)
    {
        free(I.element_nodes[i]);
    }

    free(I.C);
    free(I.node_coordinates_x);
    free(I.node_coordinates_y);
    free(I.element_nodes);
    free(I.g_H);
    free(I.g_V);
    free(I.g_num_H);
    free(I.g_num_V);
    free(I.red_nodes_H);

    fclose(I.S);
    fclose(I.R);

    for(unsigned long i = 0; i < 2*I.num_of_red_nodes_H; i++)
    {
        free(I.g_U_h[i]);
    }
    free(I.g_U_h);
    free(I.displacements);

/// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    printf("Proceeding to end it... done.\n");
    fclose(L);
    return 0;
}