#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_dv(FILE * F, double * 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%f\t]\n\n", vec[i]);
    }

    return 0;
}

int make_log_dvd(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%.12f\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_dmd(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, "%.12f\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
}

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++;
        }
    }

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

    return out_matrix;
}

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 ** vec_2_mat(double * vec, unsigned long l)
{
    double ** res = (double **)calloc(l, sizeof(double));
    for(unsigned long i = 0; i < l; i++)
    {
        res[i] = (double *)calloc(1, sizeof(double));
        res[i][0] = vec[i];
    }

    return res;
}

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

void matSub(double ** dest, double ** sour, unsigned long r, unsigned long c)
{
    for(unsigned long i = 0; i < r; i++)
    {
        for(unsigned long j = 0; j < c; j++)
        {
            dest[i][j] -= sour[i][j];
        }
    }
}

double * vecCopy(double * vec, unsigned long l)
{
    double * res = (double *)calloc(l, sizeof(double));
    for(unsigned long i = 0; i < l; i++)
    {
        res[i] = vec[i];
    }
    return res;
}

double ** matCopy(double ** mat, unsigned long r, unsigned long c)
{
    double ** res = (double **)calloc(r, sizeof(double));
    for(unsigned long i = 0; i < r; i++)
    {
        res[i] = (double *)calloc(c, sizeof(double));
        for(unsigned long j = 0; j < c; j++)
        {
            res[i][j] = mat[i][j];
        }
    }
    return res;
}

unsigned long index(unsigned long needle, unsigned long * haystack, unsigned long haystack_weight)
{
    for(unsigned long i = 0; i < haystack_weight; i++)
    {
        if(needle == haystack[i])
        {
            return i;
        }
    }
    return haystack_weight + 1;
}

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

//for(unsigned long o = 0; o < GDof; o++)
//{
//printf("sti[%lu]  %f\n", o, stiffness[o][o]);
//}

    return 1;
}

unsigned long * activate(unsigned long * contenders, unsigned long * banned, unsigned long num_con, unsigned long num_ban, unsigned long space)
{
    unsigned long num_act = num_con - num_ban;
    unsigned long * activated = (unsigned long *)calloc(2 * num_act, sizeof(unsigned long));
    unsigned long runner = 0;

    for(unsigned long i = 0; i < num_con; i++)
    {
        for(unsigned long j = 0; j < num_ban; j++)
        {
            if(contenders[i] == banned[j])
            {
                break;
            }else if(j + 1 == num_ban)
            {
                activated[runner] = i;
                activated[runner + num_act] = i + space;
                runner++;
                break;
            }
        }
    }

    return activated;
}

double ** solution(unsigned long Gdof, unsigned long * almost_active_dof, unsigned long * prescribedDof,unsigned long num_of_fixed_nodes, double ** stiffness, double * force, unsigned long num_red,/*unsigned long * &active_dof, unsigned long &num_active_dof, */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);

    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;

    active_dof = activate(almost_active_dof, prescribedDof, num_red, num_of_fixed_nodes, num_red);

    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(unsigned long i = 0; i < num_active_dof; i++)
    {
        free(buf_sti[i]);
    }
    free(buf_sti);

    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]];
            }

            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(unsigned long numberNodes, unsigned long numberElements, unsigned long ** elementNodes, double * nodeCoordinatesX, double * nodeCoordinatesY, double * displacements, double *** stresses, 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;
    FILE * S = 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 %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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]);
    }

    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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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]);
    }

    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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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 %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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] * scaleFactor), nodeCoordinatesY[i] + (displacements[i + numberNodes] * scaleFactor), z);
    }

    fprintf(S, "CELLS %lu %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;
}

bool is_in(unsigned long n, unsigned long * a, unsigned long l)
{
    for(unsigned long i = 0; i < l; i++)
    {
        if(a[i] == n)
        {
            return true;   // the number is in the array, returning index
        }
    }
    return false;   // the number is not in the array
}

unsigned long ** sub_nodes(unsigned long ** ele_nodes, unsigned long * sub_nod, unsigned long num_ele, unsigned long num_sub)
{
    unsigned long ** res = (unsigned long **)calloc(num_ele, sizeof(unsigned long));
    for(unsigned long i = 0; i < num_ele; i++)
    {
        res[i] = (unsigned long *)calloc(NODES_PER_ELEMENT, sizeof(unsigned long));
    }

    for(unsigned long i = 0; i < num_ele; i++)  // all the elements
    {
        for(unsigned long j = 0; j < NODES_PER_ELEMENT; j++)    // all the nodes
        {
            for(unsigned long x = 0; x < num_sub; x++)  // all the numberings
            {
                if(ele_nodes[i][j] == sub_nod[x])
                {
                    res[i][j] = x;
                    break;
                }
            }
        }
    }

    return res;
}

double * sub_force(unsigned long * force_indices, double ** force, unsigned long * sub_nodes, unsigned long gdof, unsigned long num_of_sub_nodes, unsigned long num_of_forces, double for_mul)
{
    double * res = (double *)calloc(gdof, sizeof(double));

    for(unsigned long i = 0; i < gdof; i++)
    {
        for(unsigned long j = 0; j < num_of_sub_nodes; j++)
        {
            for(unsigned long f = 0; f < num_of_forces; f++)
            {
                if(sub_nodes[j] == force_indices[f])
                {
                    res[j] = force[f][0] * for_mul;
                    res[j + (gdof / 2)] = force[f][1] * for_mul;
                    break;
                }
            }
        }
    }

    return res;
}

unsigned long sub_supports(unsigned long * &res, unsigned long * support_nodes, unsigned long * sub_nodes, unsigned long num_of_supps, unsigned long num_of_nodes)
{
    unsigned long supps = 0;

    for(unsigned long s = 0; s < num_of_supps; s++)
    {
        for(unsigned long n = 0; n < num_of_nodes; n++)
        {
            if(support_nodes[s] == sub_nodes[n])
            {
                supps++;
                break;
            }
        }
    }

    res = (unsigned long *)calloc(supps, sizeof(unsigned long));
    supps = 0;

    for(unsigned long s = 0; s < num_of_supps; s++)
    {
        for(unsigned long n = 0; n < num_of_nodes; n++)
        {
            if(support_nodes[s] == sub_nodes[n])
            {
                res[supps] = n;
                supps++;
                break;
            }
        }
    }

    return supps;
}

void sub_bounds(unsigned long * &sub_H, unsigned long * &sub_V, unsigned long * sub_nodes, unsigned long * H, unsigned long * V, unsigned long number_H, unsigned long number_V)
{
    sub_H = (unsigned long *)calloc(number_H, sizeof(unsigned long));
    sub_V = (unsigned long *)calloc(number_V, sizeof(unsigned long));

    for(unsigned long i = 0; i < number_H; i++) // all the H nodes
    {
        for(unsigned long j = 0; j < number_H + number_V; j++)  // all the nodes
        {
            if(H[i] == sub_nodes[j])
            {
                sub_H[i] = j;
                break;
            }
        }
    }

    for(unsigned long i = 0; i < number_V; i++) // all the V nodes
    {
        for(unsigned long j = 0; j < number_H + number_V; j++)  // all the nodes
        {
            if(V[i] == sub_nodes[j])
            {
                sub_V[i] = j;
                break;
            }
        }
    }
}

void reorder(double ** &K, double * &F, unsigned long * sub_ele_nodes, unsigned long * V, unsigned long * H, unsigned long number_V, unsigned long number_H)
{
    unsigned long nn = number_H + number_V; // number_of_nodes
    double ** k_buf = (double **)calloc(2 * nn, sizeof(double));
    double * f_buf = (double *)calloc(2 * nn, sizeof(double));
    unsigned long * HV = (unsigned long *)calloc(2 * nn, sizeof(double));

    for(unsigned long b = 0; b < nn; b++)   // creating index array later used for swapping K to submatrices
    {
        if(b < number_H)
        {
            HV[b] = H[b];
            HV[b + nn] = H[b] + nn;
        }else
        {
            HV[b] = V[b - number_H];
            HV[b + nn] = V[b - number_H] + nn;
        }

        k_buf[b] = (double *)calloc(2 * nn, sizeof(double));
        k_buf[b + nn] = (double *)calloc(2 * nn, sizeof(double));
    }

    for(unsigned long r = 0; r < 2 * nn; r++)   // copying values from K to k_buf in the swapped order
    {
        f_buf[r] = F[HV[r]];
        for(unsigned long c = 0; c < 2 * nn; c++)
        {
            k_buf[r][c] = K[HV[r]][HV[c]];
        }
    }

    for(unsigned long i = 0; i < 2 * nn; i++)
    {
        F[i] = f_buf[i];
        for(unsigned long j = 0; j < 2 * nn; j++)
        {
            K[i][j] = k_buf[i][j];
        }
    }

    for(unsigned long b = 0; b < 2 * nn; b++)
    {
        free(k_buf[b]);
    }
    free(k_buf);
    free(f_buf);
    free(HV);
}

double ** subvector(double * vec, unsigned long s, unsigned long e, unsigned long nn)    // most likely not working
{
    unsigned long l = e - s;
    double ** res = (double **)calloc(2 * l, sizeof(double));

    for(unsigned long i = 0; i < 2 * l; i++)
    {
        res[i] = (double *)calloc(1, sizeof(double));
    }

    for(unsigned long i = 0; i < l; i++)
    {
        res[i][0] = vec[i + s];
        res[i + l][0] = vec[i + s + nn];

    }

    return res;
}

double ** subvector_special(double ** vec, unsigned long * indices, unsigned long * g_indices, unsigned long l, unsigned long g_l)
{
    double ** result = (double **)calloc(2 * l, sizeof(double));
    unsigned long b;

    for(unsigned long i = 0; i < 2 * l; i++)
    {
        result[i] = (double *)calloc(1, sizeof(double));
    }

    for(unsigned long i = 0; i < l; i++)
    {
        b = index(indices[i], g_indices, g_l);
        result[i][0] = vec[b][0];
        result[i + l][0] = vec[b + g_l][0];
    }

    return result;
}

double ** submatrix(double ** mat, unsigned long r_s, unsigned long r_e, unsigned long c_s, unsigned long c_e, unsigned long nn)
{
    unsigned long R = r_e - r_s;
    unsigned long C = c_e - c_s;

    double ** res = (double **)calloc(2 * R, sizeof(double));

    for(unsigned long i = 0; i < 2 * R; i++)
    {
        res[i] = (double *)calloc(2 * C, sizeof(double));
    }

    for(unsigned long i = 0; i < R; i++)
    {
        for(unsigned long j = 0; j < C; j++)
        {
            res[i][j] = mat[i + r_s][j + c_s];
            res[i + R][j + C] = mat[i + r_s + nn][j + c_s + nn];
            res[i][j + C] = mat[i + r_s][j + c_s + nn];
            res[i + R][j] = mat[i + r_s + nn][j + c_s];
        }
    }

    return res;
}

int quasi_itoa(unsigned long num)
{
    for(int i = 1; i <= 10; i++)
    {
        if(num < pow(10, i))
        {
            return i;
        }
    }
    return -1;
}

int store_inv(double ** imat, unsigned long r, unsigned long c, char * filename)
{
    FILE * f = fopen(filename, "wb");
    if(f == NULL)
    {
        printf("store_inv() failed to open %s\n", filename);
        return -1;
    }

    fseek(f, 0, SEEK_SET);
    fwrite(&r, sizeof(unsigned long), 1, f);
    fwrite(&c, sizeof(unsigned long), 1, f);

    for(unsigned long i = 0; i < r; i++)
    {
            fwrite(imat[i], sizeof(double), c, f);
    }

    fclose(f);
    return 0;
}

int load_inv(double ** &imat, unsigned long * &r, unsigned long * &c, char * filename)
{
    FILE * f = fopen(filename, "rb");
    if(f == NULL)
    {
        printf("load_inv() failed to open %s...\n", filename);
        return -1;
    }

    fread(r, sizeof(unsigned long), 1, f);
    fread(c, sizeof(unsigned long), 1, f);

    imat = (double **)calloc(r[0], sizeof(double));

    for(unsigned long i = 0; i < r[0]; i++)
    {
        imat[i] = (double *)calloc(c[0], sizeof(double));
        fread(imat[i], sizeof(double), c[0], f);
    }

    fclose(f);
    return 0;
}

void reduce_H(double ** &K, double * &F, unsigned long num_H, unsigned long num_V, unsigned long index, const char * filename, const char * fileext, FILE * f)
{
    unsigned long nn = num_H + num_V;   // number_of_nodes
    double ** kred = NULL, ** fred = NULL;// freed
    char * name = NULL;// freed
    char * s_num = NULL;// freed

/// K = [ K_hh ; K_hv ]     U = [ U_h ]     F = [ F_h ]
///     [ K_vh ; K_vv ]         [ U_v ]         [ F_v ]

/// the formula:
/// K_hh - K_hv * K_vv^-1 * K_vh = kred

    kred = submatrix(K, 0, num_H, 0, num_H, nn);// check // freed
    double ** k_hv = submatrix(K, 0, num_H, num_H, nn, nn);   // check // freed
    double ** k_vv = submatrix(K, num_H, nn, num_H, nn, nn);  // check // freed
    double ** k_vv_inv = inv_big(k_vv,  2*num_V, f);            // check // freed
    double ** k_vh = submatrix(K, num_H, nn, 0, num_H, nn);   // check // freed
    double ** m1 = matMul(k_hv, k_vv_inv, 2*num_V, 2*num_H, 2*num_V); //check // freed
    double ** m2 = matMul(m1, k_vh, 2*num_V, 2*num_H, 2*num_H);// check // freed
    matSub(kred, m2, 2*num_H, 2*num_H);// check

    s_num = (char *)calloc(quasi_itoa(index) + 1, sizeof(char));
    sprintf(s_num, "%lu", index);
    name = (char *)calloc(BUF_LIM, sizeof(char));
    strcpy(name, filename);
    strcat(name, s_num);
    free(s_num);// freed
    strcat(name, fileext);
    store_inv(k_vv_inv, 2 * num_V, 2 * num_V, name);
    free(name);// freed

    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        free(m1[i]);// check // freed
        free(m2[i]);// check // freed
    }
    for(unsigned long i = 0; i < 2*num_V; i++)
    {
        free(k_vv[i]);// check // freed
        free(k_vh[i]);// check // freed
    }
    free(m2); // freed
    free(m1); // freed
    free(k_vv); // freed
    free(k_vh); // freed

/// the formula:
/// F_h - K_hv * K_vv^-1 * F_v = fred

    fred = subvector(F, 0, num_H, nn);// check // freed
    double ** f_v = subvector(F, num_H, nn, nn);// check // freed
    m1 = matMul(k_hv, k_vv_inv, 2*num_V, 2*num_H, 2*num_V);// check // freed
    m2 = matMul(m1, f_v, 2*num_V, 2*num_H, 1);// check // freed
    matSub(fred, m2, 2*num_H, 1);// check

    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        free(k_hv[i]);// check // freed
        free(m1[i]);// check // freed
        free(m2[i]);// check // freed
    }

    for(unsigned long i = 0; i < 2*num_V; i++)
    {
        free(k_vv_inv[i]);// check // freed
        free(f_v[i]);// check // freed
    }

    free(m2);// check // freed
    free(m1);// check // freed
    free(k_vv_inv);// check // freed
    free(k_hv);// check // freed
    free(f_v);// check // freed

/// copy kred nd fred to K and F respectively before freeing them
    for(unsigned long i = 0; i < 2 * (num_H + num_V); i++)
    {
        free(K[i]);// check
    }
    free(K);// check
    free(F);// check
    K = (double **)calloc(2*num_H, sizeof(double));// check
    F = (double *)calloc(2*num_H, sizeof(double));// check
    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        K[i] = (double *)calloc(2*num_H, sizeof(double));// check
    }

    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        F[i] = fred[i][0];// check

        for(unsigned long j = 0; j < 2*num_H; j++)
        {
            K[i][j] = kred[i][j];// check
        }
    }

    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        free(fred[i]);// check // freed
        free(kred[i]);// check // freed
    }
    free(fred);// check // freed
    free(kred);// check // freed
}

int store_reduced(FILE * f, FILE * r, double ** stiffness, double * force, double ** stiffness_reduced, double * force_reduced, unsigned long num_H, unsigned long num_V, unsigned long * H, unsigned long * V)
{
    unsigned long nn = num_H + num_V;
    unsigned long gdof = 2 * nn;
    unsigned long gH = 2 * num_H;

/// create a register
    long int subs_written = 0;
    unsigned long * written_data = (unsigned long *)calloc(1, sizeof(unsigned long));
    unsigned long data_length = 0;
    size_t r_check;

    data_length += sizeof(unsigned long) * 2;       // num_H & num_V; unsigned long; 2
    data_length += sizeof(unsigned long) * num_H;   // H; unsigned long; num_H
    data_length += sizeof(unsigned long) * num_V;   // V; unsigned long; num_V
    data_length += sizeof(double) * gdof * gdof;    // stiffness; double; gdof x gdof
    data_length += sizeof(double) * gdof;           // force; double gdof x 1
    data_length += sizeof(double) * gH * gH;        // stiffness_reduced; double; gH x gH
    data_length += sizeof(double) * gH;             // force_reduced; double; gH x 1

    //find out what position the register is already at
    fseek(r, 0, SEEK_END);
    subs_written = ftell(r) / sizeof(unsigned long);

    if(subs_written > 0)    // if this is not the first substructure to be stored, add the number of written data to current amount
    {
        fseek(r, sizeof(unsigned long) * (subs_written - 1), SEEK_SET);  //  position so that it reads the last number
        r_check = fread(written_data, sizeof(unsigned long), 1, r);
        data_length += written_data[0];
    }

    fseek(r, 0, SEEK_END);
    fwrite(&data_length, sizeof(unsigned long), 1, r);
    fseek(f, 0, SEEK_END);

/// store num_H & num_V
    fwrite(&num_H, sizeof(unsigned long), 1, f);
    fwrite(&num_V, sizeof(unsigned long), 1, f);

/// store H
    fwrite(H, sizeof(unsigned long), num_H, f);

/// store V
    fwrite(V, sizeof(unsigned long), num_V, f);

/// store stiffness
    for(unsigned long i = 0; i < gdof; i++)
    {
        fwrite(stiffness[i], sizeof(double), gdof, f);
    }

/// store force
    fwrite(force, sizeof(double), gdof, f);

/// store stiffness_reduced
    for(unsigned long i = 0; i < gH; i++)
    {
        fwrite(stiffness_reduced[i], sizeof(double), gH, f);
    }

/// store force_reduced
    fwrite(force_reduced, sizeof(double), gH, f);

    free(written_data);

    return 0;
}

int load_reduced(FILE * f, FILE * r, unsigned long n, double ** &stiffness, double * &force, double ** &stiffness_reduced, double * &force_reduced, unsigned long &num_H, unsigned long &num_V, unsigned long * &H, unsigned long * &V, int mode)
{
    unsigned long * lu_buf = NULL;
    unsigned long * data_length = (unsigned long *)calloc(2, sizeof(unsigned long));
    unsigned long nn, gdof, gH;

/// read register
    if(n == 0)
    {
        data_length[0] = 0;
    }else
    {
        fseek(r, sizeof(unsigned long) * (n - 1), SEEK_SET);
        lu_buf = (unsigned long *)calloc(1, sizeof(unsigned long));
        fread(lu_buf, sizeof(unsigned long), 1, r);
        data_length[0] = lu_buf[0];
        free(lu_buf);
    }

    fseek(r, sizeof(unsigned long) * n, SEEK_SET);
    lu_buf = (unsigned long *)calloc(1, sizeof(unsigned long));
    fread(lu_buf, sizeof(unsigned long), 1, r);
    data_length[1] = lu_buf[0];
    free(lu_buf);

/// read num_H & num_V
    fseek(f, data_length[0], SEEK_SET); // finds the start of the section in the file
    fread(&num_H, sizeof(unsigned long), 1, f);
    fread(&num_V, sizeof(unsigned long), 1, f);

    nn = num_H + num_V;
    gdof = 2 * nn;
    gH = 2 * num_H;

/// read H
    H = (unsigned long *)calloc(num_H, sizeof(unsigned long));
    fread(H, sizeof(unsigned long), num_H, f);

/// read V
    V = (unsigned long *)calloc(num_V, sizeof(unsigned long));
    fread(V, sizeof(unsigned long), num_V, f);

/// read stiffness
    if(mode == 2)   // if global stiffness is requested
    {
        stiffness = (double **)calloc(gdof, sizeof(double));
        for(unsigned long i = 0; i < gdof; i++)
        {
            stiffness[i] = (double *)calloc(gdof, sizeof(double));
            fread(stiffness[i], sizeof(double), gdof, f);
        }
    }else
    {
        fseek(f, sizeof(double) * gdof * gdof, SEEK_CUR);
    }

/// read force
    if(mode == 2)   // if global force is requested
    {
        force = (double *)calloc(gdof, sizeof(double));
        fread(force, sizeof(double), gdof, f);
    }else
    {
        fseek(f, sizeof(double) * gdof, SEEK_CUR);
    }

/// read stiffness_reduced
    if(mode == 1)   // if reduced stiffness is requested
    {
        stiffness_reduced = (double **)calloc(gH, sizeof(double));
        for(unsigned long i = 0; i < gH; i++)
        {
            stiffness_reduced[i] = (double *)calloc(gH, sizeof(double));
            fread(stiffness_reduced[i], sizeof(double), gH, f);
        }
    }else
    {
        fseek(f, sizeof(double) * gH * gH, SEEK_CUR);
    }

/// read force_reduced
    if(mode == 1)   // if reduced force is requested
    {
        force_reduced = (double *)calloc(gH, sizeof(double));
        fread(force_reduced, sizeof(double), gH, f);
    }else
    {
        fseek(f, sizeof(double) * gH, SEEK_CUR);
    }

    free(data_length);
    return 0;
}

bool check_H(unsigned long sample, std::vector<unsigned long> vec)
{
    for(unsigned long i = 0; i < vec.size(); i++)
    {
        if(vec[i] == sample)
        {
            return true;
        }
    }
    return false;
}

void globalize_H(unsigned long * &result, unsigned long &result_length, unsigned long ** H, unsigned long * num_H, unsigned long num_sub)
{
    std::vector<unsigned long> buffer;

    for(unsigned long i = 0; i < num_sub; i++)
    {
        for(unsigned long j = 0; j < num_H[i]; j++)
        {
            if(!check_H(H[i][j], buffer))
            {
                buffer.push_back(H[i][j]);
            }
        }
    }

    std::sort(buffer.begin(), buffer.end());

    result_length = (unsigned long)buffer.size();
    result = (unsigned long *)calloc(buffer.size(), sizeof(unsigned long));

    for(unsigned long i = 0; i < buffer.size(); i++)
    {
        result[i] = buffer[i];
    }

    for(unsigned long i = 0; i < result_length; i++)
    {
        buffer.erase(buffer.end()-1);
    }
}

void global_K_h(double ** &g_K, unsigned long * red_H, unsigned long num_red_H, double *** K, unsigned long ** H, unsigned long * num_H, unsigned long num_sub)
{
    g_K = (double **)calloc(2*num_red_H, sizeof(double));

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

    for(unsigned long i = 0; i < num_sub; i++)
    {
        for(unsigned long m = 0; m < num_H[i]; m++)
        {
            for(unsigned long n = 0; n < num_H[i]; n++)
            {
                g_K[index(H[i][m], red_H, num_red_H)][index(H[i][n], red_H, num_red_H)] += K[i][m][n];
                g_K[index(H[i][m], red_H, num_red_H)][index(H[i][n], red_H, num_red_H) + num_red_H] += K[i][m][n + num_H[i]];
                g_K[index(H[i][m], red_H, num_red_H) + num_red_H][index(H[i][n], red_H, num_red_H)] += K[i][m + num_H[i]][n];
                g_K[index(H[i][m], red_H, num_red_H) + num_red_H][index(H[i][n], red_H, num_red_H) + num_red_H] += K[i][m + num_H[i]][n + num_H[i]];
            }
        }
    }
}

void global_prescribed_dof(unsigned long * &g_pdof, unsigned long * pdof, unsigned long num_pdof, unsigned long * red_H, unsigned long num_red_H)
{
    g_pdof = (unsigned long *)calloc(num_pdof, sizeof(unsigned long));

    for(unsigned long i = 0; i < num_pdof; i++)
    {
        g_pdof[i] = index(pdof[i], red_H, num_red_H);
    }
}

void global_F_h(double * &g_F, unsigned long * red_H, unsigned long num_red_H, double ** F, unsigned long ** H, unsigned long * num_H, unsigned long num_sub)
{
    g_F = (double *)calloc(2*num_red_H, sizeof(double));

    for(unsigned long i = 0; i < num_sub; i++)
    {
        for(unsigned long m = 0; m < num_H[i]; m++)
        {
            g_F[index(H[i][m], red_H, num_red_H)] = F[i][m];

            g_F[index(H[i][m], red_H, num_red_H) + num_red_H] = F[i][m + num_H[i]];
        }
    }
}

void reduce_V(double * &U_v, double ** K, double * F, double ** U_h, /*unsigned long * active, */unsigned long num_H, unsigned long num_V, /*unsigned long num_active,*/ unsigned long * H, unsigned long * g_H, unsigned long num_g_H, unsigned long index, const char * filename, const char * fileext, FILE * f)
{
    unsigned long nn = num_H + num_V;   // number_of_nodes
    char * name = NULL; // freed
    unsigned long * r = (unsigned long *)calloc(1, sizeof(unsigned long)); // freed
    unsigned long * c = (unsigned long *)calloc(1, sizeof(unsigned long)); // freed
    double ** buf_U_h = NULL; // freed

    U_v = (double *)calloc(2 * num_V, sizeof(double));
    buf_U_h = subvector_special(U_h, H, /*active, */g_H, num_H/*, num_active*/, num_g_H);

/// K = [ K_hh ; K_hv ]     U = [ U_h ]     F = [ F_h ]
///     [ K_vh ; K_vv ]         [ U_v ]         [ F_v ]

/// the formula:
/// U_v = K_v^-1 * (F_v - K_vh * U_h)

    double ** k_vh = submatrix(K, num_H, nn, 0, num_H, nn); // freed

    name = (char *)calloc(BUF_LIM - 1, sizeof(char));
    sprintf(name, "%s%lu%s", filename, index, fileext);
    double ** k_vv_inv = NULL; // freed
    load_inv(k_vv_inv, r, c, name);
    free(name); // freed

    double ** m1 = matMul(k_vh, buf_U_h, 2*num_H, 2*num_V, 1); // freed
    double ** f_v = subvector(F, num_H, nn, nn); // freed

    matSub(f_v, m1, 2*num_V, 1);
    double ** m2 = matMul(k_vv_inv, f_v, 2*num_V, 2*num_V, 1); // freed

    for(unsigned long i = 0; i < 2 * num_V; i++)
    {
        U_v[i] = m2[i][0];
    }

    for(unsigned long i = 0; i < 2*num_V; i++)
    {
        free(k_vv_inv[i]); // freed
        free(k_vh[i]); // freed
        free(m1[i]); // freed
        free(f_v); // freed
        free(m2[i]); // freed
    }
    free(m1); // freed
    free(k_vv_inv); // freed
    free(f_v); // freed
    free(k_vh); // freed
    free(m2); // freed

    for(unsigned long i = 0; i < 2*num_H; i++)
    {
        free(buf_U_h[i]); // freed
    }
    free(buf_U_h); // freed

    free(r); // freed
    free(c); // freed
}

void store_disp_H(double ** U_H, double * U, unsigned long * H, unsigned long num_H, unsigned long num_nodes)
{
    for(unsigned long i = 0; i < num_H; i++)
    {
        U[H[i]] = U_H[i][0];
        U[H[i] + num_nodes] = U_H[i + num_H][0];
    }
}

void store_disp_V(double * U_V, double * U, unsigned long * V, unsigned long num_V, unsigned long num_nodes)
{
    for(unsigned long i = 0; i < num_V; i++)
    {
        U[V[i]] = U_V[i];
        U[V[i] + num_nodes] = U_V[i + num_V];
    }
}

void load_stuff(double * &x, double * &y, unsigned long ** &el_nodes, FILE * f, unsigned long num_nodes, unsigned long num_els)
{
    char * buf = NULL;
    unsigned long node_check = 0;
    unsigned long el_check = 0;
    buf = (char *)calloc(BUF_LIM + 1, sizeof(char));

    while(!feof(f))
    {
        fgets(buf, BUF_LIM, f);

        if(strncmp(buf, ".<number_of_nodes>", 18) == 0)
        {
            node_check = strtoul(parse_input(buf, '{', '}'), NULL, 10);

            x = (double *)calloc(node_check, sizeof(double));
            y = (double *)calloc(node_check, sizeof(double));

            for(unsigned long n = 0; n < node_check; n++)
            {
                fgets(buf, BUF_LIM, f);
                x[n] = strtod(parse_input(buf, '{', ';'), NULL);
                y[n] = strtod(parse_input(buf, ';', '}'), NULL);
            }

            while(true)
            {
                fgets(buf, BUF_LIM, f);
                if(strncmp(buf, "#<number_of_elements>", 21) == 0)
                {
                    el_check = strtoul(parse_input(buf, '{', '}'), NULL, 10);
                    el_nodes = (unsigned long **)calloc(el_check, sizeof(unsigned long));
                    for(unsigned long e = 0; e < el_check; e++)
                    {
                        el_nodes[e] = (unsigned long *)calloc(NODES_PER_ELEMENT, sizeof(unsigned long));

                        fgets(buf, BUF_LIM, f);
                        el_nodes[e][0] = strtoul(parse_input(buf, '{', ','), NULL, 10);
                        el_nodes[e][1] = strtoul(parse_input(buf, ',', ':'), NULL, 10);
                        el_nodes[e][2] = strtoul(parse_input(buf, ':', ';'), NULL, 10);
                        el_nodes[e][3] = strtoul(parse_input(buf, ';', '}'), NULL, 10);
                    }
                    break;
                }
            }
            break;
        }
    }
}

#endif // HEAD_H_INCLUDED
