#ifndef HEAD_H_INCLUDED
#define HEAD_H_INCLUDED

//#define IN_FUNCTION_LOGS 1
//#define IN_FUNCTION_LOGS_DETAILED 2

int setup_gen() //freed
{
    FILE * s;

    s = fopen("setup.txt", "w");
    if(s != NULL)
    {
        char msg[] = "<E>={220000000000}\n<nu>={0.3}\n<force_multiplier>={1000}\n<density>={7850}\n<thickness>={0.001}\n<scale_factor>={100}";  // the default settings to be written in the newly generated file
        fwrite(msg, sizeof(msg) - 1, 1, s); // writing the message
        fclose(s);  // closing the file
    }
    free(s);
    return 0;
}

char equals[] = " = ";

int make_log_lu(FILE * F, unsigned long num, char msg[])    //nothing to be freed
{
    fprintf(F, "%s%s%lu\n\n", msg, equals, num);

    return 0;
}

int make_log_lu(FILE * F, unsigned long num, char msg[], unsigned long index)   //nothing to be freed
{
    fprintf(F, "%s[%lu]%s%lu\n\n", msg, index, equals, num);

    return 0;
}

int make_log_d(FILE * F, double num, char msg[])    //nothing to be freed
{
    fprintf(F, "%s%s%f\n\n", msg, equals, num);

    return 0;
}

int make_log_d(FILE * F, double num, char msg[], unsigned long index)   //nothing to be freed
{
    fprintf(F, "%s[%lu]%s%f\n\n", msg, index, equals, num);

    return 0;
}

int make_log_coor(FILE * F, double * x, double * y, unsigned long r, char msg[])    //nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[\t%f\t%f\t]\n\n", x[i], y[i]);
    }

    return 0;
}

int make_log_dv(FILE * F, double * vec, unsigned long c, char msg[])    //nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);
    for(unsigned long i = 0; i < c; i++)
    {
        fprintf(F, "[\t%f\t]\n\n", vec[i]);
    }

    return 0;
}

int make_log_luv(FILE * F, unsigned long * vec, unsigned long c, char msg[])    //nothing to be freed
{

    fprintf(F, "%s%s\n", msg, equals);
    for(unsigned long i = 0; i < c; i++)
    {
        fprintf(F, "[\t%lu\t]\n\n", vec[i]);
    }

    return 0;
}

int make_log_luv(FILE * F, unsigned long * vec, unsigned long c, char msg[], unsigned long index)   //nothing to be freed
{

    fprintf(F, "%s[%lu]%s\n", msg, index, equals);
    for(unsigned long i = 0; i < c; i++)
    {
        fprintf(F, "[\t%lu\t]\n\n", vec[i]);
    }

    return 0;
}

int make_log_lum(FILE * F, unsigned long ** mat, unsigned long r, unsigned long c, char msg[])  //nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[\t");
        for(unsigned long j = 0; j < c; j++)
        {
            fprintf(F, "%lu\t", mat[i][j]);
        }
        fprintf(F, "]\n\n");
    }

    return 0;
}

int make_log_dm(FILE * F, double ** mat, unsigned long r, unsigned long c, char msg[])  //nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[\t");
        for(unsigned long j = 0; j < c; j++)
        {
            fprintf(F, "%f\t", mat[i][j]);
        }
        fprintf(F, "]\n\n");
    }

    return 0;
}

int make_log_dm(FILE * F, double ** mat, unsigned long r, unsigned long c, char msg[], unsigned long index) //nothing to be freed
{
    fprintf(F, "%s[%lu]%s\n", msg, index, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[\t");
        for(unsigned long j = 0; j < c; j++)
        {
            fprintf(F, "%f\t", mat[i][j]);
        }
        fprintf(F, "]\n\n");
    }

    return 0;
}

int make_log_dmf(FILE * F, double ** mat, unsigned long r, unsigned long c, char msg[])  //nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[\t");
        for(unsigned long j = 0; j < c; j++)
        {
            fprintf(F, "%.9f\t", mat[i][j]);
        }
        fprintf(F, "]\n\n");
    }

    return 0;
}

int make_log_dt(FILE * F, double *** ten, unsigned long r, unsigned long c, unsigned long d, char msg[])    // nothing to be freed
{
    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < r; i++)
    {
        fprintf(F, "[ ");
        for(unsigned long j = 0; j < c; j++)
        {
            fprintf(F, "{ ");
            for(unsigned long k = 0; k < d; k++)
            {
                fprintf(F, " %.9f ", ten[i][j][k]);
            }
            fprintf(F, " } ");
        }
        fprintf(F, "]\n\n");
    }

    return 0;
}

int make_log_reactions(FILE * F, double ** rea, unsigned long * prescribed_dof, unsigned long num_of_prescribed_dofs, unsigned long num_of_nodes, char msg[])
{
    double x_sum = 0, y_sum = 0;

    fprintf(F, "%s%s\n", msg, equals);

    for(unsigned long i = 0; i < num_of_prescribed_dofs; i++)
    {
        fprintf(F, "%lu_x:\t[ %f ]\n\n", prescribed_dof[i], rea[prescribed_dof[i]][0]);
        x_sum += rea[prescribed_dof[i]][0];
    }
    for(unsigned long i = 0; i < num_of_prescribed_dofs; i++)
    {
        fprintf(F, "%lu_y:\t[ %f ]\n\n", prescribed_dof[i], rea[prescribed_dof[i] + num_of_nodes][0]);
        y_sum += rea[prescribed_dof[i] + num_of_nodes][0];
    }

    fprintf(F, "X_sum = %f;\nY_sum = %f\n\n", x_sum, y_sum);

    return 0;
}

double ** makeC(double E, double nu)    // generates C matrix for stress strain calculation // is freed afterwards
{
    double ** C;
    C = (double **)calloc(3, sizeof(double));   // 3 because of 3x3 matrix

    for(unsigned short i = 0; i < 3; i++)
    {
        C[i] = (double *)calloc(3, sizeof(double)); //
    }

    C[0][0] = E / (1 - pow(nu,2));
    C[0][1] = (nu * E) / (1 - pow(nu,2));
    C[0][2] = 0;

    C[1][0] = C[0][1];
    C[1][1] = C[0][0];
    C[1][2] = 0;

    C[2][0] = 0;
    C[2][1] = 0;
    C[2][2] = E / (2 * (1 + nu));

    return C;
}

char * parse_input(char * buf, char d1, char d2)   // given a string and delimeters, returns the name of the variable   //nothing to be freed
{
    short s = -1;    // indicator for the start of variable name
    short e = -1;    // indicator for the end of variable name
    char * output = NULL;   // buffer for the name to be returned

    for(unsigned long i = 0; i < strlen(buf); i++)  // goes through the whole line
    {
        if(buf[i] == d1)    // position of start of the variable name
        {
            s = i;
        }else if(buf[i] == d2)  // position of start of the variable name
        {
            e = i;
        }

        if((s != -1) && (e != -1))  // if both limits have been found
        {
            break;
        }
    }

    if((s == -1) || (e == -1))
    {
        return NULL;
    }

    output = (char*)calloc((e - s), sizeof(char));  // allocate space needed for the output buffer

    for(short i = 0; i < e - s - 1; i++)    // loop for going through the bounds of the parsed string
    {
        output[i] = buf[i + s + 1]; // extracting only the correct characters
    }
    output[e - s] = '\0';   // terminating the string

    return output;  // returning the string for further use
}

void run_memory_run(variables &I)   // clears all the memory before terminating the program
{
    unsigned long i, j;
    if(I.stresses != NULL)
    {
        for(i = 0; i < I.number_of_elements; i++)
        {
            for(j = 0; j < NODES_PER_ELEMENT; j++)
            {
                free(I.stresses[i][j]);
            }
            free(I.stresses[i]);
        }
        free(I.stresses);
    }
    if(I.support_nodes != NULL)
    {
        //printf("support_nodes\n");
        for(i = 0; i < I.number_of_support_types; i++)
        {
            free(I.support_nodes[i]);
        //    printf("%lu ", i);
        }
        free(I.support_nodes);
        //printf("\n");
    }
    if(I.support_names != 0)
    {
        //printf("support_names\n");
        for(i = 0; i < I.number_of_support_types; i++)
        {
            free(I.support_names[i]);
        //    printf("%lu ", i);
        }
        free(I.support_names);
        //printf("\n");
    }
    if(I.element_nodes != NULL)
    {
        //printf("element_nodes\n");
        for(i = 0; i < I.number_of_elements; i++)
        {
            free(I.element_nodes[i]);
        //    printf("%lu ", i);
        }
        free(I.element_nodes);
        //printf("\n");
    }
    if(I.substructure_elements != NULL)
    {
        //printf("substructure_elements\n");
        for(i = 0; i < I.number_of_substructures; i++)
        {
            free(I.substructure_elements[i]);
        //    printf("%lu ", i);
        }
        free(I.substructure_elements);
        //printf("\n");
    }
    if(I.forces != NULL)
    {
        //printf("forces\n");
        for(i = 0; i < I.number_of_forces; i++)
        {
            free(I.forces[i]);
        //    printf("%lu ", i);
        }
        free(I.forces);
        //printf("\n");
    }
    if(I.stiffness != NULL)
    {
        //printf("stiffness\n");
        for(i = 0; i < I.number_of_nodes; i++)
        {
            free(I.stiffness[i]);
            //printf("%lu ", i);
        }
        free(I.stiffness);    // need to finish
        //printf("\n");
    }
    if(I.displacements != NULL)
    {
        //printf("displacements\n");
        for(i = 0; i < I.global_dof; i++)
        {
            free(I.displacements[i]);
            //printf("%lu ", i);
        }
        free(I.displacements);
        //printf("\n");
    }
    if(I.mass != NULL)
    {
        //printf("mass\n");
        for(i = 0; i < I.number_of_nodes; i++)
        {
            free(I.mass[i]);
            //printf("%lu ", i);
        }
        free(I.mass);
        //printf("\n");
    }
    if(I.weights != NULL)
    {
        //printf("weights\n");
        for(i = 0; i < 4; i++)
        {
            free(I.weights[i]);
            //printf("%lu ", i);
        }
        free(I.weights);
        //printf("\n");
    }
    if(I.locations != NULL)
    {
        //printf("locations\n");
        for(i = 0; i < 4; i++)
        {
            free(I.locations[i]);
            //printf("%lu ", i);
        }
        free(I.locations);
        //printf("\n");
    }
    if(I.C != NULL)
    {
        //printf("C\n");
        for(i = 0; i < 3; i++)
        {
            free(I.C[i]);
            //printf("%lu ", i);
        }
        free(I.C);
        //printf("\n");
    }
    if(I.fil_buf != NULL)
    {
        free(I.fil_buf);
        //printf("fil_buf\n");
    }
    if(I.name_buf != NULL)
    {
        free(I.name_buf);
        //printf("name_buf\n");
    }
    if(I.var_buf != NULL)
    {
        free(I.var_buf);
        //printf("var_buf\n");
    }
    if(I.support_index)
    {
        free(I.support_index);
        //printf("support_index\n");
    }
    if(I.number_of_support_nodes_per_type != NULL)
    {
        free(I.number_of_support_nodes_per_type);
        //printf("number_of_support_nodes_per_type\n");
    }
    if(I.number_of_substructure_elements != NULL)
    {
        free(I.number_of_substructure_elements);
        //printf("number_of_substructure_elements\n");
    }
    if(I.force_indices != NULL)
    {
        free(I.force_indices);
        //printf("force_indices\n");
    }
    if(I.node_coordinates_x != NULL)
    {
        free(I.node_coordinates_x);
        //printf("node_coordinates_x\n");
    }
    if(I.node_coordinates_y != NULL)
    {
        free(I.node_coordinates_y);
        //printf("node_coordinates_y\n");
    }
    if(I.force != NULL)
    {
        free(I.force);
        //printf("force\n");
    }
    if(I.reactions != NULL)
    {
        free(I.reactions);
        //printf("reactions\n");
    }
    if(I.F != NULL)
    {
        free(I.F);
        //printf("F\n");
    }
    if(I.S != NULL)
    {
        free(I.S);
    }
}

double ** gaussQuadratureLoc()  //nothing to be freed
{
    double ** locations = NULL;
    locations = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));

    for(unsigned short i = 0; i < NODES_PER_ELEMENT; i++)
    {
        locations[i] = (double *)calloc(2, sizeof(double));
    }

    locations[0][0] = -0.577350269189626;
    locations[0][1] = -0.577350269189626;

    locations[1][0] = 0.577350269189626;
    locations[1][1] = -0.577350269189626;

    locations[2][0] = 0.577350269189626;
    locations[2][1] = 0.577350269189626;

    locations[3][0] = -0.577350269189626;
    locations[3][1] = 0.577350269189626;

    return locations;
}

double ** gaussQuadratureWei()  //nothing to be freed
{
    double ** weights = NULL;
    weights = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));

    for(unsigned short i = 0; i < NODES_PER_ELEMENT; i++)
    {
        weights[i] = (double *)calloc(1, sizeof(double));
    }

    weights[0][0] = 1;
    weights[1][0] = 1;
    weights[2][0] = 1;
    weights[3][0] = 1;

    return weights;
}

double ** shapeFunctionQ4shape(double xi, double eta)   //nothing to be freed
{
    double ** shape = NULL;

    shape = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));

    for(unsigned short i = 0; i < NODES_PER_ELEMENT; i++)
    {
        shape[i] = (double *)calloc(1, sizeof(double));
    }

    shape[0][0] = ((1 - xi) * (1 - eta)) / 4;
    shape[1][0] = ((1 + xi) * (1 - eta)) / 4;
    shape[2][0] = ((1 + xi) * (1 + eta)) / 4;
    shape[3][0] = ((1 - xi) * (1 + eta)) / 4;

    return shape;
}

double ** shapeFunctionQ4naturalDerivatives(double xi, double eta)  //nothing to be freed
{
    double ** derivatives = NULL;

    derivatives = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));

    for(unsigned short i = 0; i < NODES_PER_ELEMENT; i++)
    {
        derivatives[i] = (double *)calloc(2, sizeof(double));
    }

    derivatives[0][0] = -(1 - eta) / 4;
    derivatives[0][1] = -(1 - xi) / 4;

    derivatives[1][0] = (1 - eta) / 4;
    derivatives[1][1] = -(1 + xi) / 4;

    derivatives[2][0] = (1 + eta) / 4;
    derivatives[2][1] = (1 + xi) / 4;

    derivatives[3][0] = -(1 + eta) / 4;
    derivatives[3][1] = (1 - xi) / 4;

    return derivatives;
}

double determinant(double **mat, unsigned long k, int print_status = 1)  // freed
{
    double s = 1, det = 0;
    double ** buf;
    unsigned long i,j,m,n,c;

    buf = (double **)calloc(k, sizeof(double));
    for(i = 0; i < k; i++)
    {
        buf[i] = (double *)calloc(k, sizeof(double));
    }

    if(k == 1)
    {
        return (mat[0][0]);
    }
    else
    {
        det = 0;

        for(c = 0; c < k; c++)
        {
            m = 0;
            n = 0;

            for(i = 0; i < k; i++)
            {
                for(j = 0; j < k; j++)
                {
                    buf[i][j] = 0;
                    if(i != 0 && j != c)
                    {
                        buf[m][n] = mat[i][j];
                        if(n < (k-2))
                        {
                            n++;
                        }
                        else
                        {
                            n=0;
                            m++;
                        }
                    }
                }
            }

            det = det + s * (mat[0][c] * determinant(buf, k - 1, 0));
            s = -1 * s;
        }
    }

    for(i = 0; i < k; i++)
    {
        free(buf[i]);
    }
    free(buf);

    return det;
}

double ** transpose(double ** num, double ** fac, unsigned long r)  //freed
{
    unsigned long i, j;
    double ** b, ** inverse, d;

    inverse = (double **)calloc(r, sizeof(double));
    b = (double **)calloc(r, sizeof(double));
    for(i = 0; i < r; i++)
    {
        inverse[i] = (double *)calloc(r, sizeof(double));
        b[i] = (double *)calloc(r, sizeof(double));
    }

    for(i = 0; i < r; i++)
    {
        for(j = 0; j < r; j++)
        {
            b[i][j] = fac[j][i];
        }
    }

    d = determinant(num,r);

    for(i = 0; i < r; i++)
    {
        for(j = 0; j < r; j++)
        {
            inverse[i][j] = b[i][j] / d;
        }
    }

    for(i = 0; i < r; i++)
    {
        free(b[i]);
    }
    free(b);

    return inverse;
}

double ** cofactor(double ** num, unsigned long f)  //freed
{
    double ** b, ** fac;
    unsigned long p,q,m,n,i,j;

    b = (double **)calloc(f, sizeof(double));
    fac = (double **)calloc(f, sizeof(double));
    for(i = 0; i < f; i++)
    {
        b[i] = (double *)calloc(f, sizeof(double));
        fac[i] = (double *)calloc(f, sizeof(double));
    }

    for(q = 0; q < f; q++)
    {
        for(p = 0; p < f; p++)
        {
            m = 0;
            n = 0;
            for(i = 0; i < f; i++)
            {
                for(j = 0; j < f; j++)
                {
                    if(i != q && j != p)
                    {
                        b[m][n] = num[i][j];
                        if(n < (f - 2))
                        {
                            n++;
                        }
                        else
                        {
                            n=0;
                            m++;
                        }
                    }
                }
            }

        fac[q][p] = pow(-1, q + p) * determinant(b, f - 1, 0);
        }
    }

    for(i = 0; i < f; i++)
    {
        free(b[i]);
    }
    free(b);

    return transpose(num, fac, f);
}

double ** inv(double ** a, unsigned long k, int print_status = 1) //nothing to be freed
{
    double d;

    d = determinant(a, k, print_status);

    if(d == 0)
    {
        printf("\nInverse of Entered Matrix is not possible\n");
    }
    else
    {
        return cofactor(a, k);
    }

    return NULL;
}

double ** inv_big(double ** in_matrix, unsigned long n, FILE * F)
{
    double rat;
    double a;
    unsigned long l;
    double ** matrix = NULL;
    double ** out_matrix = NULL;

    matrix = (double **)calloc(n, sizeof(double));
    out_matrix = (double **)calloc(n, sizeof(double));

    for(unsigned long i = 0; i < n; i++)
    {
        matrix[i] = (double *)calloc(2 * n, sizeof(double));   // 2n because the algorithm requires it
        out_matrix[i] = (double *)calloc(n, sizeof(double));

        for(unsigned long j = 0; j < n; j++)
        {
            matrix[i][j] = in_matrix[i][j]; // copying the input matrix to the one that will be used in the algorithm
        }
    }

    for(unsigned long i = 0; i < n; i++)
    {
        for(unsigned long j = n; j < 2 * n; j++)
        {
            if(i == (j - n))
            {
                matrix[i][j] = 1;
            }else
            {
                matrix[i][j] = 0;
            }
        }
    }

    for(unsigned long i = 0; i < n; i++)
    {
        for(unsigned long j = 0; j < n; j++)
        {
            if(i != j)
            {
                rat = matrix[j][i] / matrix[i][i];

                for(unsigned long k = 0; k < 2 * n; k++)
                {
                    matrix[j][k] -= rat * matrix[i][k];
                }
            }
        }
    }

    for(unsigned long i = 0; i < n; i++)
    {
        a = matrix[i][i];

        for(unsigned long j = 0; j < 2 * n; j++)
        {
            matrix[i][j] /= a;
        }
    }

    for(unsigned long i = 0; i < n; i++)
    {
        l = 0;

        for(unsigned long j = n; j < 2 * n; j++)
        {
            out_matrix[i][l] = matrix[i][j];
            l++;
        }
    }

    return out_matrix;
}

double ** LU(double ** mat, unsigned long n)    // the non working lu
{

    unsigned long i, j, k;
	double x;

    double ** M = NULL;
    M = (double **)calloc(n, sizeof(double));
    for(i = 0; i < n; i++)
    {
        M[i] = (double *)calloc(n, sizeof(double));

        for(j = 0; j < n; j++)
        {
            M[i][j] = mat[i][j];    // copy the matrix to perform decomposition on
        }
    }

    for(k = 0; k <= n - 1; k++)
    {
        for(j = k + 1; j <= n; j++)
        {
            x = M[j][k] / M[k][k];

            for(i = k; i <= n; i++)
            {
                M[j][i] = M[j][i] - x * M[k][i];
            }
            M[j][k] = x;
        }
    }

    return M;
}

double ** inv(double ** M, unsigned long n, FILE * L)   // the non working lu inverse
{
    unsigned long i, j, m;
    double x;

    double ** lu = NULL;
    double ** s = NULL;
    double * d = NULL;
    double * y = NULL;

    s = (double **)calloc(n, sizeof(double));
    d = (double *)calloc(n, sizeof(double));
    y = (double *)calloc(n, sizeof(double));

    for(i = 0; i < n; i++)
    {
        s[i] = (double *)calloc(n, sizeof(double));
    }

    lu = LU(M, n);

    for(m=0;m<=2;m++)
    {
        for(i = 0; i < n; i++)
        {
            d[i] = 0;
        }
        d[m]=1.0;

        for(i=0;i<=n;i++)
        {
            x=0;

            for(j=0;j<=i-1;j++)
            {
                x=x+lu[i][j]*y[j];
            }

            y[i]=(d[i]-x);
        }

        for(i=n;i>=0;i--)
        {
            x=0;

            for(j=i+1;j<=n;j++)
            {
                x=x+lu[i][j]*s[j][m];
            }

            s[i][m]=(y[i]-x)/lu[i][i];
        }
	}

	make_log_dm(L, s, n, n, "inverse s");

	for(i = 0; i < n; i++)
    {
        free(lu[i]);
    }
    free(lu);
    free(d);
    free(y);

    return s;
}

double ** jacobianMatrix(double ** elementCoordinates, double ** derivatives)//nothing to be freed
{
    double ** jac = NULL;

    jac = (double **)calloc(2, sizeof(double));
    for(unsigned short i = 0; i < 2; i++)
    {
        jac[i] = (double *)calloc(2, sizeof(double));
    }

    for(unsigned short n = 0; n < NODES_PER_ELEMENT; n++)   // matrix multiplication with transposed node coordinates
    {
        for(unsigned short i = 0; i < 2; i++)
        {
            for(unsigned short j = 0; j < 2; j++)
            {
                jac[i][j] += elementCoordinates[n][i] * derivatives[n][j];
            }
        }
    }

    return jac;
}

double ** XYD(double ** derivatives, double ** invJac)  //nothing to be freed
{
    double ** xyd = NULL;

    xyd = (double **)calloc(4, sizeof(double));
    for(unsigned short i = 0; i < 4; i++)
    {
        xyd[i] = (double *)calloc(2, sizeof(double));
    }

    for(unsigned short n = 0; n < 2; n++)
    {
        for(unsigned short i = 0; i < 4; i++)
        {
            for(unsigned short j = 0; j < 2; j++)
            {
                xyd[i][j] += derivatives[i][n] * invJac[n][j];
            }
        }
    }
    return xyd;
}

double jacobianDet(double ** jacob) //nothing to be freed
{
    double determinant;

    determinant = (jacob[0][0] * jacob[1][1]) - (jacob[0][1] * jacob[1][0]);

    return determinant;
}

double ** matMul(double ** m1, double ** m2, unsigned long M, unsigned long r1, unsigned long c2, int print_status = 1)   //nothing to be freed
{
    double ** result = NULL;
    result = (double **)calloc(r1, sizeof(double));
    for(unsigned long q = 0; q < r1; q++)
    {
        result[q] = (double *)calloc(c2, sizeof(double));
    }

    for(unsigned long m = 0; m < M; m++)
    {
        if(print_status == 0)
        {
            printf("progress: %lu / %lu\n", m, M);
        }

        for(unsigned long i = 0; i < r1; i++)
        {
            for(unsigned long j = 0; j < c2; j++)
            {
                result[i][j] += m1[i][m] * m2[m][j];
            }
        }
    }

    return result;
}

double ** matMul_firstTransposed(double ** m1, double ** m2, unsigned long R, unsigned long c1, unsigned long c2)   //nothing to be freed
{
    double ** result = NULL;
    result = (double **)calloc(c1, sizeof(double));
    for(unsigned long q = 0; q < c1; q++)
    {
        result[q] = (double *)calloc(c2, sizeof(double));
    }

    for(unsigned long m = 0; m < R; m++)
    {
        for(unsigned long i = 0; i < c1; i++)
        {
            for(unsigned long j = 0; j < c2; j++)
            {
                result[i][j] += m1[m][i] * m2[m][j];
            }
        }
    }

    return result;
}

double ** matMul_scalar(double ** m, double s, unsigned long r, unsigned long c)    //nothing to be freed
{
    double ** result = NULL;
    result = (double **)calloc(r, sizeof(double));
    for(unsigned long i = 0; i < r; i++)
    {
        result[i] = (double *)calloc(c, sizeof(double));
    }

    for(unsigned long i = 0; i < r; i++)
    {
        for(unsigned long j = 0; j < c; j++)
        {
            result[i][j] = m[i][j] * s;
        }
    }

    return result;
}

int formStiffness2D(double ** &stiffness, double ** &mass, unsigned long GDof, unsigned long numberElements, unsigned long ** elementNodes, unsigned long numberNodes, double * nodeCoordinatesX, double * nodeCoordinatesY, double ** C, double rho, double thickness, double ** gaussLocations, double ** gaussWeights, FILE * F)
{
    unsigned long indice[4];
    unsigned long elementDof[8];
    unsigned long i, j, e, q;
    double ** element_coordinates = NULL;   //freed
    double xi = 0;
    double eta = 0;
    double ** shape = NULL; //freed
    double ** derivatives = NULL;   //freed
    double ** jacob = NULL; //freed
    double ** invJacob = NULL;  //freed
    double ** XYderivatives = NULL; //freed
    double ** B_mat = NULL; //freed
    double detJac = 0;
    double ** buf_matrix1 = NULL;
    double ** buf_matrix2 = NULL;
    double ** buf_matrix3 = NULL;
    double ** buf_matrix4 = NULL;
    double ** buf_matrix5 = NULL;

    stiffness = (double **)calloc(GDof, sizeof(double));
    mass = (double **)calloc(GDof, sizeof(double));
    for(i = 0; i < GDof; i++)
    {
        stiffness[i] = (double *)calloc(GDof, sizeof(double));
        mass[i] = (double *)calloc(GDof, sizeof(double));
    }

    element_coordinates = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));
    for(i = 0; i < NODES_PER_ELEMENT; i++)
    {
        element_coordinates[i] = (double *)calloc(2, sizeof(double));
    }

    for(e = 0; e < numberElements; e++)
    {
        for(i = 0; i < NODES_PER_ELEMENT; i++)
        {
            indice[i] = elementNodes[e][i];
        }

        for(i = 0; i < NODES_PER_ELEMENT; i++)
        {
            elementDof[i] = indice[i];
            elementDof[i + NODES_PER_ELEMENT] = indice[i] + numberNodes;
        }

        for(q = 0; q < NODES_PER_ELEMENT; q++)
        {
            xi = gaussLocations[q][0];
            eta = gaussLocations[q][1];

            shape = shapeFunctionQ4shape(xi, eta);
            derivatives = shapeFunctionQ4naturalDerivatives(xi, eta);

            for(i = 0; i < 4; i++)
            {
                free(shape[i]);
            }
            free(shape);

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                element_coordinates[i][0] = nodeCoordinatesX[indice[i]];
                element_coordinates[i][1] = nodeCoordinatesY[indice[i]];
            }

            jacob = jacobianMatrix(element_coordinates, derivatives);
            invJacob = inv(jacob, 2/*, F*/);
            XYderivatives = XYD(derivatives, invJacob);

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                free(derivatives[i]);
            }
            free(derivatives);

            for(i = 0; i < 2; i++)
            {
                free(invJacob[i]);
            }
            free(invJacob);

            detJac = jacobianDet(jacob);

            for(i = 0; i < 2; i++)
            {
                free(jacob[i]);
            }
            free(jacob);

            B_mat = (double **)calloc(3, sizeof(double));
            for(i = 0; i < 3; i++)
            {
                B_mat[i] = (double *)calloc(2 * NODES_PER_ELEMENT, sizeof(double));
            }

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                B_mat[0][i] = XYderivatives[i][0];
                B_mat[1][i + NODES_PER_ELEMENT] = XYderivatives[i][1];
                B_mat[2][i] = XYderivatives[i][1];
                B_mat[2][i + NODES_PER_ELEMENT] = XYderivatives[i][0];
            }

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                free(XYderivatives[i]);
            }
            free(XYderivatives);

            buf_matrix1 = matMul_firstTransposed(B_mat, C, 3, 8, 3);
            buf_matrix2 = matMul_scalar(buf_matrix1, thickness, 8, 3);

            for(i = 0; i < 8; i++)
            {
                free(buf_matrix1[i]);
            }
            free(buf_matrix1);

            buf_matrix3 = matMul(buf_matrix2, B_mat, 3, 8, 8);

            for(i = 0; i < 8; i++)
            {
                free(buf_matrix2[i]);
            }
            free(buf_matrix2);

            buf_matrix4 = matMul_scalar(buf_matrix3, gaussWeights[q][0], 8, 8);

            for(i = 0; i < 3; i++)
            {
                free(B_mat[i]);
            }
            free(B_mat);


            for(i = 0; i < 8; i++)
            {
                free(buf_matrix3[i]);
            }
            free(buf_matrix3);

            buf_matrix5 = matMul_scalar(buf_matrix4, detJac, 8, 8);

            for(i = 0; i < 8; i++)
            {
                free(buf_matrix4[i]);
            }
            free(buf_matrix4);

            for(i = 0; i < 2 * NODES_PER_ELEMENT; i++)
            {
                for(j = 0; j < 2 * NODES_PER_ELEMENT; j++)
                {
                    stiffness[elementDof[i]][elementDof[j]] += buf_matrix5[i][j];
                }
            }

            for(i = 0; i < 8; i++)
            {
                free(buf_matrix5[i]);
            }
            free(buf_matrix5);
        }
    }

    for(i = 0; i < NODES_PER_ELEMENT; i++)
    {
        free(element_coordinates[i]);
    }
    free(element_coordinates);

    return 1;
}

double ** solution(unsigned long Gdof, unsigned long * prescribedDof,unsigned long num_of_fixed_nodes, double ** stiffness, double * force, FILE * F)   // freed
{
    double ** dis = NULL;   // cannot be freed
    double ** buf_sti = NULL;   // freed
    double ** buf_force = NULL; // freed
    double ** buf_dis = NULL;   // freed
    unsigned long * active_dof = NULL;  // freed
    unsigned long i = 0, j = 0;
    unsigned long searched_one = 0;

    unsigned long num_active_dof = Gdof - (2 * num_of_fixed_nodes);

    active_dof = (unsigned long *)calloc(num_active_dof, sizeof(unsigned long));

    buf_sti = (double **)calloc(num_active_dof, sizeof(double));
    buf_force = (double **)calloc(num_active_dof, sizeof(double));
    for(unsigned long i = 0; i < num_active_dof; i++)
    {
        buf_sti[i] = (double *)calloc(num_active_dof, sizeof(double));
        buf_force[i] = (double *)calloc(1, sizeof(double));
    }

    i = 0;
    j = 0;

    while(i < num_active_dof)   // generation of active_dof array
    {
        j = 0;

        while(j < num_of_fixed_nodes)
        {
            if(searched_one == prescribedDof[j])
            {
                searched_one++;
                break;
            }else if(searched_one == (prescribedDof[j] + (Gdof / 2)))
            {
                searched_one++;
                break;
            }else if((j + 1) == num_of_fixed_nodes)
            {
                active_dof[i] = searched_one;
                i++;
                searched_one++;
            }
            j++;
        }
    }

    for(i = 0; i < num_active_dof; i++)
    {
        for(j = 0; j < num_active_dof; j++)
        {
            buf_sti[i][j] = stiffness[active_dof[i]][active_dof[j]];
        }
        buf_force[i][0] = force[active_dof[i]];
    }

    buf_dis = matMul(inv_big(buf_sti, num_active_dof, F), buf_force, num_active_dof, num_active_dof, 1);

    for(i = 0; i < num_active_dof; i++)
    {
        free(buf_force[i]);
    }
    free(buf_force);

    dis = (double **)calloc(Gdof, sizeof(double));
    for(i = 0; i < Gdof; i++)
    {
        dis[i] = (double *)calloc(1, sizeof(double));
    }

    for(i = 0; i < num_active_dof; i++)
    {
            dis[active_dof[i]][0] = buf_dis[i][0];
    }

    for(i = 0; i < num_active_dof; i++)
    {
        free(buf_dis[i]);
    }
    free(buf_dis);
    free(active_dof);

    return dis;
}

double *** stresses2D(unsigned long Gdof, unsigned long numberElements, unsigned long ** elementNodes, unsigned long numberNodes, double * nodeCoordinatesX, double * nodeCoordinatesY, double ** displacements, double ** C)
{
    unsigned long i, j, e, q;
    double ** gaussWeights = NULL;  // 4x1, freed
    double ** gaussLocations = NULL;    // 4x2, freed
    double *** stress = NULL;   // numberElements x NODES_PER_ELEMENT x 3, cannot be freed
    double ** elementCoordinates = NULL;    // 4x2, freed
    unsigned long indice[NODES_PER_ELEMENT];
    unsigned long elementDof[2 * NODES_PER_ELEMENT];
    double ** stress_buf;   // 3x1, freed
    double ** dis_buf;  // 8x1, freed
    double xi = 0;
    double eta = 0;
    double ** naturalDerivatives = NULL;    //4x2, freed
    double ** jac = NULL;   // 2x2, freed
    double ** invJac = NULL;    // 2x2, freed
    double ** XYDerivatives = NULL; // 4x2, freed
    double ** B_mat = NULL; //3 x 2*NODES_PER_ELEMENT,freed

    gaussWeights = gaussQuadratureWei();
    gaussLocations = gaussQuadratureLoc();

    stress = (double ***)calloc(numberElements, sizeof(double));
    for(i = 0; i < numberElements; i++)
    {
        stress[i] = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));
        for(j = 0; j < NODES_PER_ELEMENT; j++)
        {
            stress[i][j] =(double *)calloc(3, sizeof(double));
        }
    }

    elementCoordinates = (double **)calloc(NODES_PER_ELEMENT, sizeof(double));
    for(i = 0; i < NODES_PER_ELEMENT; i++)
    {
        elementCoordinates[i] = (double *)calloc(2, sizeof(double));
    }

    stress_buf = (double **)calloc(3, sizeof(double));
    for(i = 0; i < 3; i++)
    {
        stress_buf[i] = (double *)calloc(1, sizeof(double));
    }

    dis_buf = (double **)calloc(2 * NODES_PER_ELEMENT, sizeof(double));
    for(i = 0; i < 2 * NODES_PER_ELEMENT; i++)
    {
        dis_buf[i] = (double *)calloc(1, sizeof(double));
    }

    B_mat = (double **)calloc(3, sizeof(double));
    for(i = 0; i < 3; i++)
    {
        B_mat[i] = (double *)calloc(2 * NODES_PER_ELEMENT, sizeof(double));
    }

    for(e = 0; e < numberElements; e++)
    {
        for(i = 0; i < NODES_PER_ELEMENT; i++)
        {
            indice[i] = elementNodes[e][i];
            elementDof[i] = indice[i];
            elementDof[i + NODES_PER_ELEMENT] = indice[i] + numberNodes;
        }

        for(q = 0; q < NODES_PER_ELEMENT; q++)
        {
            xi = gaussLocations[q][0];
            eta = gaussLocations[q][1];

            naturalDerivatives = shapeFunctionQ4naturalDerivatives(xi, eta);

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                elementCoordinates[i][0] = nodeCoordinatesX[indice[i]];
                elementCoordinates[i][1] = nodeCoordinatesY[indice[i]];
            }

            jac = jacobianMatrix(elementCoordinates, naturalDerivatives);
            invJac = inv(jac, 2);
            XYDerivatives = XYD(naturalDerivatives, invJac);

            for(i = 0; i < NODES_PER_ELEMENT; i++)
            {
                B_mat[0][i] = XYDerivatives[i][0];
                B_mat[1][i + NODES_PER_ELEMENT] = XYDerivatives[i][1];
                B_mat[2][i] = XYDerivatives[i][1];
                B_mat[2][i + NODES_PER_ELEMENT] = XYDerivatives[i][0];
            }

            for(i = 0; i < 2 * NODES_PER_ELEMENT; i++)
            {
                dis_buf[i][0] = displacements[elementDof[i]][0];
            }

            stress_buf = matMul(C, matMul(B_mat, dis_buf, 8, 3, 1), 3, 3, 1);

            for(i = 0; i < 3; i++)
            {
                stress[e][q][i] = stress_buf[i][0];
            }
        }
    }

    for(i = 0; i < 2; i++)
    {
        free(jac[i]);
        free(invJac[i]);
    }
    free(jac);
    free(invJac);

    for(i = 0; i < 3; i++)
    {
        free(stress_buf[i]);
        free(B_mat[i]);
    }
    free(stress_buf);
    free(B_mat);

    for(i = 0; i < NODES_PER_ELEMENT; i++)
    {
        free(gaussWeights[i]);
        free(gaussLocations[i]);
        free(elementCoordinates[i]);

        free(dis_buf[2 * i]);
        free(dis_buf[2 * i + 1]);

        free(naturalDerivatives[i]);
        free(XYDerivatives[i]);
    }
    free(gaussWeights);
    free(gaussLocations);
    free(elementCoordinates);
    free(dis_buf);
    free(naturalDerivatives);
    free(XYDerivatives);

    return stress;
}

int make_solutions(FILE * S, unsigned long numberNodes, unsigned long numberElements, unsigned long ** elementNodes, double * nodeCoordinatesX, double * nodeCoordinatesY, double ** displacements, double *** stresses, double ** reactions, double scaleFactor, char timestamp[])
{
    double ** stress_sum = NULL;
    double sigma1, sigma2, centre, radius;
    unsigned long i, j, k;
    double z = 0.0;
    unsigned char * num_of_memberships = NULL;

    stress_sum = (double **)calloc(numberNodes, sizeof(double));
    num_of_memberships = (unsigned char *)calloc(numberNodes, sizeof(unsigned char));
    for(i = 0; i < numberNodes; i++)
    {
        stress_sum[i] = (double *)calloc(3, sizeof(double));
    }

    for(i = 0; i < numberElements; i++)
    {
        for(j = 0; j < NODES_PER_ELEMENT; j++)
        {
            for(k = 0; k < 3; k++)
            {
                stress_sum[elementNodes[i][j]][k] += stresses[i][j][k];
            }

            num_of_memberships[elementNodes[i][j]] += 1;
        }
    }

    for(i = 0; i < numberNodes; i++)
    {
        for(j = 0; j < 3; j++)
        {
            stress_sum[i][j] = stress_sum[i][j] / num_of_memberships[i];
        }
    }

    S = fopen("solutions/solution/undeformed.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/undeformed.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "undeformed state of the model, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%f %f %f\n", nodeCoordinatesX[i], nodeCoordinatesY[i], z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fclose(S);

    S = fopen("solutions/solution/deformed.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_displacements_x.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_displacements_x.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with x component of displacements, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS displacement_x double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", displacements[i][0]);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_displacements_y.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_displacements_y.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with y component of displacements, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS displacement_y double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", displacements[i + numberNodes][0]);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_stress_sigma_x.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_stresses_sigma_x.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with x component of stress, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS sigma_x double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", stress_sum[i][0]);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_stress_sigma_y.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_stresses_sigma_y.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with y component of stress, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS sigma_y double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", stress_sum[i][1]);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_stress_tau_xy.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_stresses_tau_xy.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with xy component of shear stress, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS tau_xy double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", stress_sum[i][2]);
    }

    fclose(S);

    S = fopen("solutions/solution/undeformed_indices.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/undeformed_indices.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "undeformed state of the model with node indices, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i], nodeCoordinatesY[i], z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS undeformed_indices unsigned_long\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%lu\n", i);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_indices.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_indices.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with node indices, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS deformed_indices unsigned_long\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%lu\n", i);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_stress_tresca.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_stresses_tresca.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with Tresca stress, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS sigma_Tresca double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        centre = (stress_sum[i][0] + stress_sum[i][1]) / 2;
        radius = sqrt(pow(((stress_sum[i][0] - stress_sum[i][1]) / 2), 2) + pow(stress_sum[i][2], 2));

        sigma1 = centre + radius;
        sigma2 = centre - radius;

        if((sigma1 > 0) || (sigma2 < 0))
        {
            fprintf(S, "%.9f\n", sigma1 - sigma2);
        }else if((sigma1 < 0) || (sigma2 < 0))
        {
            fprintf(S, "%.9f\n", -sigma2);
        }else
        {
            fprintf(S, "%.9f\n", sigma1);
        }
    }

    fclose(S);

    for(i = 0; i < numberNodes; i++)
    {
        free(stress_sum[i]);
    }
    free(stress_sum);

    S = fopen("solutions/solution/deformed_reactions_x.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_reactions_x.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with x component of reaction forces, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS reactions_x double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", reactions[i][0]);
    }

    fclose(S);

    S = fopen("solutions/solution/deformed_reactions_y.vtu", "w");

    if(S == NULL)
    {
        printf("Cannot open solutions/solution/deformed_reactions_y.vtu... Make sure the folder exists and repeat the simulation.\n");
        return -1;
    }

    fprintf(S, "# vtk DataFile Version 3. 0.\n");
    fprintf(S, "deformed state of the model with y component of reaction forces, generated: %s", timestamp);
    fprintf(S, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
    fprintf(S, "POINTS %lu double\n", numberNodes);

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f %.9f %.9f\n", nodeCoordinatesX[i] + (displacements[i][0] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes][0] * scaleFactor), z);
    }

    fprintf(S, "CELLS %d %lu\n", numberElements, numberElements * (NODES_PER_ELEMENT + 1));
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "%d %lu %lu %lu %lu\n", NODES_PER_ELEMENT, elementNodes[i][0], elementNodes[i][1], elementNodes[i][2], elementNodes[i][3]);
    }

    fprintf(S, "CELL_TYPES %lu\n", numberElements);
    for(i = 0; i < numberElements; i++)
    {
        fprintf(S, "9\n");
    }

    fprintf(S, "POINT_DATA %lu\n", numberNodes);
    fprintf(S, "SCALARS reactions_y double\n");
    fprintf(S, "LOOKUP_TABLE default\n");

    for(i = 0; i < numberNodes; i++)
    {
        fprintf(S, "%.9f\n", reactions[i + numberNodes][0]);
    }

    fclose(S);

    return 0;
}

#endif // HEAD_H_INCLUDED
