/*
    This program performs finite element analysis
    by reading in the data, doing assembling, and
    then solving the linear system for a truss

	Updated 7/4/99

    SLFFEA source file
    Version:  1.0
    Copyright (C) 1999  San Le 

    The source code contained in this file is released under the
    terms of the GNU Library General Public License.
 
*/

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "tsconst.h"
#include "tsstruct.h"

int tswriter ( BOUND , int *, double *, int *, double *, int *, MATL *,
        char *, STRAIN *, STRESS *, double *);

int solve(double *,double *,int *, int );

int decomp(double *, int *, int );

int tsKassemble(double *, int *, double *, int *, double *, int *, int *,
	int *, MATL *, STRAIN *, STRESS *, double *);

int diag( int *, int *, int, int, int, int);

int formlm( int *, int *, int *, int, int, int );

int tsformid( BOUND, int *);

int tsreader( BOUND , int *, double *, int *, double *, MATL *, FILE *,
	STRESS *, double *);

int dof, analysis_flag, neqn, nmat, numel, numnp, sof;

int main(int argc, char** argv)
{
	int i,j;
        int *id, *lm, *idiag, check;
        XYZI *mem_XYZI;
        int *mem_int, sofmi, sofmf, sofmXYZI, ptr_inc;
        MATL *matl;
        double *mem_double;
        double fpointx, fpointy, fpointz;
	int *connect, *el_matl, dum;
	double *coord, *force, *U, *A;
	char name[20];
	char buf[ BUFSIZ ];
        FILE *o1;
	BOUND bc;
        STRESS *stress;
        STRAIN *strain;
        long timec;
        long timef;

        sof = sizeof(double);

	bzero(name,20*sizeof(char));

    	printf("What is the name of the file containing the \n");
    	printf("truss structural data? \n");
    	scanf( "%20s",&name);

/*   o1 contains all the structural data */

        o1 = fopen( name ,"r" );

	if(o1 == NULL) {
		printf("error on open\n");
		exit(1);
	}
	fgets( buf, BUFSIZ, o1 );
        fscanf( o1, "%d %d %d\n ",&numel,&numnp,&nmat);
        dof=numnp*ndof;

/*   Begin allocation of meomory */

/* For the doubles */
        sofmf=2*dof+numnp*nsd;
        mem_double=(double *)calloc(sofmf,sizeof(double));

        if(!mem_double )
        {
                printf( "failed to allocate memory for double\n ");
                exit(1);
        }
                                           ptr_inc = 0;
        coord=(mem_double+ptr_inc);        ptr_inc += numnp*nsd;
        force=(mem_double+ptr_inc);        ptr_inc += dof;
        U=(mem_double+ptr_inc);            ptr_inc += dof;

/* For the materials */
        matl=(MATL *)calloc(nmat,sizeof(MATL));
        if(!matl )
        {
                printf( "failed to allocate memory for matl doubles\n ");
                exit(1);
        }

/* For the STRESS doubles */

        stress=(STRESS *)calloc(numel,sizeof(STRESS));
        if(!stress )
        {
                printf( "failed to allocate memory for stress doubles\n ");
                exit(1);
        }

/* For the STRAIN doubles */

        strain=(STRAIN *)calloc(numel,sizeof(STRAIN));
        if(!strain )
        {
                printf( "failed to allocate memory for strain doubles\n ");
                exit(1);
        }

/* For the integers */
        sofmi= numel*npel+2*dof+numel*npel*ndof+numel+numnp+1;
        mem_int=(int *)calloc(sofmi,sizeof(int));
        if(!mem_int )
        {
                printf( "failed to allocate memory for integers\n ");
                exit(1);
        }
                                                ptr_inc = 0;
        connect=(mem_int+ptr_inc);              ptr_inc += numel*npel;
        id=(mem_int+ptr_inc);                   ptr_inc += dof;
        idiag=(mem_int+ptr_inc);                ptr_inc += dof;
        lm=(mem_int+ptr_inc);                   ptr_inc += numel*npel*ndof;
        el_matl=(mem_int+ptr_inc);              ptr_inc += numel;
	bc.force =(mem_int+ptr_inc);            ptr_inc += numnp;
	bc.num_force=(mem_int+ptr_inc);         ptr_inc += 1;

/* For the XYZI integers */
        sofmXYZI=numnp+1;
        mem_XYZI=(XYZI *)calloc(sofmXYZI,sizeof(XYZI));
        if(!mem_XYZI )
        {
                printf( "failed to allocate memory for XYZI integers\n ");
                exit(1);
        }
                                             ptr_inc = 0;
        bc.fix =(mem_XYZI+ptr_inc);          ptr_inc += numnp;
        bc.num_fix=(mem_XYZI+ptr_inc);       ptr_inc += 1;

        timec = clock();
        timef = 0;

	check = tsreader( bc, connect, coord, el_matl, force, matl, o1, stress, U);
        if(!check) printf( " Problems with tsreader \n");

        check = tsformid( bc, id );
        if(!check) printf( " Problems with tsformid \n");

#if DATA_ON
        printf( "\n This is the id matrix \n");
        for( i = 0; i < numnp; ++i )
        {
                printf("\n node(%4d)",i);
                for( j = 0; j < ndof; ++j )
                {
                        printf(" %4d  ",*(id+ndof*i+j));
                }
        }
#endif

	check = formlm( connect, id, lm, ndof, npel, numel );
        if(!check) printf( " Problems with formlm \n");

#if DATA_ON
        printf( "\n\n This is the lm matrix \n");
        for( i = 0; i < numel; ++i )
        {
            printf("\n element(%4d)",i);
            for( j = 0; j < neqel; ++j )
            {
                printf( "%5d   ",*(lm+neqel*i+j));
            }
        }
        printf( "\n");
#endif

        check = diag( idiag, lm, ndof, neqn, npel, numel);
        if(!check) printf( " Problems with diag \n");

#if DATA_ON
        printf( "\n\n This is the idiag matrix \n");
        for( j = 0; j < neqn; ++j )
        {
            printf( "\ndof %5d   %5d",j,*(idiag+j));
        }
        printf( "\n");
#endif

/*   allocate meomory for A, the skyline representation of the global stiffness */

        A=(double *)calloc((*(idiag+neqn-1)+1),sizeof(double));

	analysis_flag = 1;
        check = tsKassemble( A, connect, coord, el_matl, force, id, idiag, lm,
		matl, strain, stress, U);
        if(!check) printf( " Problems with assembler \n");

#if DATA_ON
        printf( "\n\n This is the force matrix \n");
        for( j = 0; j < neqn; ++j )
        {
            printf( "\ndof %5d   %14.5f",j,*(force+j));
        }
        printf(" \n");
#endif

/* Perform LU Crout decompostion on the system */

        check = decomp(A,idiag,neqn);
        if(!check) printf( " Problems with decomp \n");

/* Solve the system */

        check = solve(A,force,idiag,neqn);
        if(!check) printf( " Problems with solve \n");
#if DATA_ON
        printf( "\n This is the solution to the problem \n");
#endif
        for( i = 0; i < dof; ++i )
        {
                if( *(id + i) > -1 )
                {
                        *(U + i) = *(force + *(id + i));
                }
        }
#if DATA_ON
        for( i = 0; i < numnp; ++i )
        {
                if( *(id+ndof*i) > -1 )
                {
                	printf("\n node %3d x   %14.6f ",i,*(U+ndof*i));
                }
                if( *(id+ndof*i+1) > -1 )
                {
                	printf("\n node %3d y   %14.6f ",i,*(U+ndof*i+1));
                }
                if( *(id+ndof*i+2) > -1 )
                {
                	printf("\n node %3d z   %14.6f ",i,*(U+ndof*i+2));
                }
        }
        printf(" \n");
#endif

/* Calculate the reaction forces */
	analysis_flag = 2;
	bzero(force,dof*sof);
#if DATA_ON
        printf( "\n\n These are the axial displacements and forces \n");
#endif
        check = tsKassemble( A, connect, coord, el_matl, force, id, idiag, lm,
		matl, strain, stress, U);
        if(!check) printf( " Problems with assembler \n");
#if DATA_ON
        printf( "\n\n These are the reaction forces \n");
        for( i = 0; i < numnp; ++i )
        {
                if( *(id+ndof*i) < 0 )
                {
                	printf("\n node %3d x   %14.6f ",i,*(force+ndof*i));
                }
                if( *(id+ndof*i+1) < 0 )
                {
                	printf("\n node %3d y   %14.6f ",i,*(force+ndof*i+1));
                }
                if( *(id+ndof*i+2) < 0 )
                {
                	printf("\n node %3d z   %14.6f ",i,*(force+ndof*i+2));
                }
        }

        printf( "\n\n               These are the updated coordinates \n");
        printf( "\n                  x               y             z \n");

        for( i = 0; i < numnp; ++i )
        {
                fpointx = *(coord+nsd*i) + *(U+ndof*i);
                fpointy = *(coord+nsd*i+1) + *(U+ndof*i+1);
                fpointz = *(coord+nsd*i+2) + *(U+ndof*i+2);
                printf("\n node %3d %14.9f %14.9f %14.9f",i,fpointx,fpointy,fpointz);
        }
	printf(" \n");
#endif
	check = tswriter ( bc, connect, coord, el_matl, force, id, matl,
	        name, strain, stress, U);
	if(!check) printf( " Problems with writer \n");

	timec = clock();
        printf("\n elapsed CPU = %lf\n\n",( (double)timec)/800.);

	free(mem_double);
	free(strain);
	free(stress);
	free(mem_int);
	free(mem_XYZI);
}