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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

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

        FILE * S = NULL;

        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 ** displacements = NULL;
        double ** mass = NULL;
        double * node_coordinates_x = NULL;
        double * node_coordinates_y = NULL;
        double * force = NULL;
        double ** reactions = NULL;

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

        double *** stresses = NULL;

        unsigned long fixed_index = 0;

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

        unsigned long num_buf = 0;  // used to buffer numerical values for comparison purposes
    };

#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
        run_memory_run(I);  // free the memory
        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
        run_memory_run(I);  // free the memory
        return -1;
    }

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

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

    I.F = fopen("model.mkp", "r");    // opening the model.mkp file with geometry info
    if(I.F == NULL)   // if it did not open
    {
        printf("model.mkp not found.\nMake sure the file is located in the same folder as this program's executable.\n");    // error message
        run_memory_run(I);  // free the memory
        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

                I.node_coordinates_x = (double *)calloc(I.number_of_nodes, sizeof(double));  // allocate memory for x coordinates
                I.node_coordinates_y = (double *)calloc(I.number_of_nodes, sizeof(double));  // allocate memory for y coordinates
                I.global_dof = 2 * I.number_of_nodes;   // set the number of global degrees of freedom

                for(unsigned long i = 0; i < I.number_of_nodes; i++)  // goes through all the nodes
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F); // reads the line

                    I.node_coordinates_x[i] = strtod(parse_input(I.fil_buf, '{', ';'), NULL);   // and stores the x coordinate
                    I.node_coordinates_y[i] = strtod(parse_input(I.fil_buf, ';', '}'), NULL);   // together with the y coordinate
                }

            }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

                I.element_nodes = (unsigned long **)calloc(I.number_of_elements, sizeof(unsigned long));    // allocate memory for storing node indices of all the elements

                for(unsigned long i = 0; i < I.number_of_elements; i++)   // goes through all the elements
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F); // read the line

                    I.element_nodes[i] = (unsigned long *)calloc(NODES_PER_ELEMENT, sizeof(unsigned long));   // allocate memory for node indices of a given element

                    I.element_nodes[i][0] = strtoul(parse_input(I.fil_buf, '{', ','), NULL, 10);    // store the node indices
                    I.element_nodes[i][1] = strtoul(parse_input(I.fil_buf, ',', ':'), NULL, 10);    // hard coded, because why not??
                    I.element_nodes[i][2] = strtoul(parse_input(I.fil_buf, ':', ';'), NULL, 10);
                    I.element_nodes[i][3] = strtoul(parse_input(I.fil_buf, ';', '}'), NULL, 10);
                }
            }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);

                I.number_of_substructure_elements = (unsigned long *)calloc(I.number_of_substructures, sizeof(unsigned long));  // allocates memory for storing amounts of elements in substructures
                I.substructure_elements = (unsigned long **)calloc(I.number_of_substructures, sizeof(unsigned long));   // allocate memory for storing indices of elements in substructures

                for(unsigned long i = 0; i < I.number_of_substructures; i++)   // goes through all the substructures
                {
                    fgets(I.fil_buf, BUF_LIM - 1, I.F);

                    I.number_of_substructure_elements[i] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);

                    I.substructure_elements[i] = (unsigned long *)calloc(I.number_of_substructure_elements[i], sizeof(unsigned long));  // gets the amount of elements in the substructure and allocates memory needed

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

                        I.substructure_elements[i][j] = strtoul(parse_input(I.fil_buf, '{', '}'), NULL, 10);    // store the index of element in the substructure
                    }
                }
            }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
                                run_memory_run(I);  // free the memory
                                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(I.number_of_forces, 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
    }

    I.locations = gaussQuadratureLoc();
    I.weights = gaussQuadratureWei();

printf("about to make stiffness\n");
    formStiffness2D(I.stiffness, I.mass, I.global_dof, I.number_of_elements, I.element_nodes, I.number_of_nodes, I.node_coordinates_x, I.node_coordinates_y, I.C, I.density, I.thickness, I.locations, I.weights, L);
printf("just made stiffness\n");

    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.force = (double *)calloc(I.global_dof, sizeof(double));   // allocating force vectors
    for(unsigned long i = 0; i < I.number_of_forces; i++)
    {
        I.force[I.force_indices[i]] =  I.force_multiplier * I.forces[i][0];    // the  component
        I.force[I.force_indices[i] + I.number_of_nodes] =  I.force_multiplier * I.forces[i][1];    // the y component
    }

printf("about to make displacements\n");
    I.displacements = solution(I.global_dof, I.support_nodes[I.fixed_index], I.number_of_support_nodes_per_type[I.fixed_index], I.stiffness, I.force, L);
    I.reactions = matMul(I.stiffness, I.displacements, I.global_dof, I.global_dof, 1);
printf("just made displacements\n");

    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);

printf("Simulation done.\n");
    fclose(L);

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

    run_memory_run(I);  // free the memory
    return 0;
}