#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>
#include "types.h"
#include "mesh.h"

/*-------------------------------------------------------------------------*/
/* Function prototypes                                                     */
int
getLine( FILE *fp, 	/* File pointer for input                    */
         char *s, 	/* Character buffer to be loaded             */
         int  lim  	/* Buffer size in characters                 */);

void
readHeader(FILE *fp);

void
readCoord(FILE *fpi, FILE *fpo);

void
readConn(FILE *fpi, FILE *fpo);

void
ReadMeshNodeset( FILE *fpi, FILE *fpo );

void
ReadMeshSideset( FILE *fpi, FILE *fpo );

void
echoHeader(FILE *fp);

int
iscomment( char *cdat );

void
Allocate();

int
bandwidth();

void
formknf( const Real T,     /* Wire tension */
	 const int ltype,  /* load type =0 constant, =1 split */
	 Real *k,          /* stiffness */
	 Real *f );        /* rhs (force */

void
blastbc( Real *k, Real *f );

void
write_u( Real *u, FILE *fpo, FILE *fpg);

int
udu( Real *a,        /* coefficient matrix in udu format */
     Real *y,        /* RHS on input, solution on output */
     const int neq,  /* No. of equations to solve        */
     const int ib,   /* 1/2 bandwidth of the system      */
     const int lu    /* lu=1 factorize, lu=0 resolve     */ );

/* Function prototypes                                                     */
/*-------------------------------------------------------------------------*/

/* Global data */
MeshData mesh;
Oned2    *oned;     /* These are allocatable objects */
Quad4    *quad4;
Hex8     *hex;
Nodeset  *ndset;    
Sideset  *sdset;    
/* Global data */

int
main(int argc, char **argv)
/***********************************************************************

Program: fem1d - 1D FEM harness

Author : Mark A. Christon

Date   : 01/06/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  FILE *fpi,*fpo,*fpg;
  Real *k, *f;
  Real T,w0;
  int  i, nbw, ltype;
 
  if (argc<4) {
   fprintf(stdout,"\n\tUsage: fem1d mesh output plot\n\n");
   exit(0);
  } 

  /* Open input & output files */
  fpi=fopen(argv[1],"r");
  fpo=fopen(argv[2],"w");
  fpg=fopen(argv[3],"w");

  /* Read & echo header data */
  readHeader(fpi);
  echoHeader(fpo);

  /* Allocate memory */
  Allocate();

  /* Read & Echo nodal coordinates, connectivity and nodesets */
  readCoord(fpi,fpo);
  readConn(fpi,fpo);
  if (mesh.Nnd_sets>0) ReadMeshNodeset(fpi,fpo);

  /* Compute the 1/2-bandwidth & allocate for coefficient matrix, rhs */
  nbw=bandwidth();
  k=(Real *)malloc(mesh.Nnp*nbw*sizeof(Real));
  f=(Real *)malloc(mesh.Nnp*   sizeof(Real));

  /* Form K and F */
  bzero(k,sizeof(Real)*mesh.Nnp*nbw);
  bzero(f,sizeof(Real)*mesh.Nnp   );
  T=30000.0;
  ltype=0;
  fprintf(stdout,"\t1/2-bandwidth= %4d\n",nbw);
  fprintf(stdout,"\tTension      = %10.4e\n",T);
  fprintf(stdout,"\tLoad type    = %4d\n",ltype);
  fprintf(fpo   ,"\n");
  fprintf(fpo   ,"\t1/2-bandwidth= %4d\n",nbw);
  fprintf(fpo   ,"\tTension      = %10.4e\n",T);
  fprintf(fpo   ,"\tLoad type    = %4d\n",ltype);
  formknf(T,ltype,k,f);
 
  /* Set EBC's based on nodesets */
  blastbc(k,f);

  /* Solve Ku=F */
  if (udu(k,f,mesh.Nnp,nbw,1)<0) {
    fprintf(stdout,"\n\tError termindation from linear solver ...\n");
    exit(0);
  }

  /* Output displacements */
  write_u(f,fpo,fpg);

  /* Close all files */
  fclose(fpi);
  fclose(fpo);
  fclose(fpg);
}

void
Allocate()
/***********************************************************************

Routine: Allocate - allocate for coordinates, connectivity, nodesets

Author : Mark A. Christon

Date   : 01/07/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{

  switch (mesh.Ndim) {
    case ONED:
      oned=(Oned2 *)malloc(mesh.Nel2*sizeof(Oned2));
      mesh.x[0]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      break;

    case TWOD:
      quad4 =(Quad4 *)malloc(mesh.Nel4*sizeof(Quad4));
      mesh.x[0]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      mesh.x[1]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      break;

    case THREED:
      hex =(Hex8 *)malloc(mesh.Nel8*sizeof(Hex8));
      mesh.x[0]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      mesh.x[1]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      mesh.x[2]=(Real *)malloc(mesh.Nnp*sizeof(Real));
      break;

    default:
      fprintf(stdout,"Allocate: Problem in element type ...\n");
      exit(0);
      break;
  }

  /* Allocate nodesets */
  if(mesh.Nnd_sets>0) ndset=(Nodeset *)malloc(mesh.Nnd_sets*sizeof(Nodeset));
  if(mesh.Nsd_sets>0) sdset=(Sideset *)malloc(mesh.Nsd_sets*sizeof(Sideset));
}

void
readHeader(FILE *fp)
/***********************************************************************

Routine: readHeader - read the mesh header information

Author : Mark A. Christon

Date   : 01/06/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int Nch,iset;
  int Nnp,Nel,Nnpe,Nmat;
  char cbuf[MAXCHR];
  char ctok[10];
  mesh.Ndim=Nel=Nmat=0;
  mesh.Nnd_sets=mesh.Nsd_sets=0;

  /* Read the title */
  if (getLine(fp,mesh.title,MAXCHR)<=0) {
    fprintf(stdout,"readHeader: Problem reading title line...\n");
    exit(0);
  }

  /* Parse the header */
  iset=0;
  while (iset==0) {
    Nch=getLine(fp,cbuf,80);
    if (strncasecmp("#",cbuf,1)==0) {
      iset=0;
    }
    else if (strncasecmp("Nnd_sets",cbuf,8)==0) {
      sscanf(cbuf,"%s %d",ctok,&mesh.Nnd_sets);
    }
    else if (strncasecmp("Nsd_sets",cbuf,8)==0) {
      sscanf(cbuf,"%s %d",ctok,&mesh.Nsd_sets);
    }
    else if (strncasecmp("Ndim",cbuf,4)==0) {
      sscanf(cbuf,"%s %d",ctok,&mesh.Ndim);
    }
    else if (strncasecmp("Nmat",cbuf,4)==0) {
      sscanf(cbuf,"%s %d",ctok,&Nmat);
    }
    else if (strncasecmp("Nnpe",cbuf,4)==0) {
      sscanf(cbuf,"%s %d",ctok,&Nnpe);
    }
    else if (strncasecmp("Nnp",cbuf,3)==0) {
      sscanf(cbuf,"%s %d",ctok,&mesh.Nnp);
    }
    else if (strncasecmp("Nel",cbuf,3)==0) {
      sscanf(cbuf,"%s %d",ctok,&Nel);
    }
    else if (strncasecmp("end",cbuf,3)==0) {
      iset=1;
    }
    else {
      fprintf(stdout,"Unknown command ...\n");
      fprintf(stdout,"%s\n",cbuf);
    }
  }

  /* Setup parameters according to dimension */
  switch (mesh.Ndim) {
    case ONED:
      	mesh.Nel2 =Nel;
      	mesh.Nnpe2=Nnpe;
        mesh.Nmat2=Nmat;
	break;

    case TWOD:
      	mesh.Nel4 =Nel;
      	mesh.Nnpe4=Nnpe;
        mesh.Nmat4=Nmat;
	break;

    case THREED:
      	mesh.Nel8 =Nel;
      	mesh.Nnpe8=Nnpe;
        mesh.Nmat8=Nmat;
	break;

    default:
	fprintf(stdout,"readHead: Problem in element type ...\n");
 	exit(0);
	break;
  }
}

void
readConn(FILE *fpi, FILE *fpo)
/***********************************************************************

Routine: readConn - read connectivity

Author : Mark A. Christon

Date   : 01/06/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int iset;
  int i,mat,n1,n2,n3,n4,n5,n6,n7,n8;
  char cbuf[MAXCHR];
  char ctok[10];

  /* Parse connectivity */
  fprintf(fpo,"\n\tElement Connectivity    ");
  fprintf(fpo,"\n\t====================\n\n");
  switch (mesh.Ndim) {

  case ONED:
    fprintf(fpo,"\tElement node-1 node-2\n");
    fprintf(fpo,"\t------- ------ ------\n");
    for (i=0;i<mesh.Nel2;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%5d%10d%10d",&mat,&n1,&n2);
          oned[i].mat=mat;
          oned[i].ix[0]=n1;
          oned[i].ix[1]=n2; 
	  fprintf(fpo,"\t%7d %6d %6d\n",
		  oned[i].mat,oned[i].ix[0],oned[i].ix[1]);
	  iset=1;
	}
#ifdef DEBUG
	else {
	  fprintf(stdout,"%s\n",cbuf);
	}
#endif
      }
    }
    break;

  case TWOD:
    fprintf(fpo,"\tElement node-1 node-2 node-3 node-4\n");
    fprintf(fpo,"\t------- ------ ------ ------ ------\n");
    for (i=0;i<mesh.Nel4;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%5d%10d%10d%10d%10d",&mat,&n1,&n2,&n3,&n4);
          quad4[i].mat=mat;
          quad4[i].ix[0]=n1;
          quad4[i].ix[1]=n2; 
          quad4[i].ix[2]=n3; 
          quad4[i].ix[3]=n4; 
	  fprintf(fpo,"\t%7d %6d %6d %6d %6d\n",
		  quad4[i].mat,quad4[i].ix[0],quad4[i].ix[1],
		               quad4[i].ix[2],quad4[i].ix[3]);
	  iset=1;
	}
      }
    }
    break;

  case THREED:
    fprintf(fpo,"\tElement node-1 node-2 node-3 node-4");
    fprintf(fpo,"node-5 node-6 node-7 node-8");
    fprintf(fpo,"\t------- ------ ------ ------ ------");
    fprintf(fpo,"------ ------ ------ ------\n");
    for (i=0;i<mesh.Nel8;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%5d%10d%10d%10d%10d%10d%10d%10d%10d",
                 &mat,&n1,&n2,&n3,&n4,&n5,&n6,&n7,&n8);
          hex[i].mat=mat;
          hex[i].ix[0]=n1;
          hex[i].ix[1]=n2; 
          hex[i].ix[2]=n3; 
          hex[i].ix[3]=n4; 
          hex[i].ix[4]=n5; 
          hex[i].ix[5]=n6; 
          hex[i].ix[6]=n7; 
          hex[i].ix[7]=n8; 
	  fprintf(fpo,"\t%7d %6d %6d %6d %6d %6d %6d %6d %6d\n",
		  hex[i].mat,
		  hex[i].ix[0],hex[i].ix[1],hex[i].ix[2],hex[i].ix[3],
		  hex[i].ix[4],hex[i].ix[5],hex[i].ix[6],hex[i].ix[7]);
	  iset=1;
	}
      }
    }
    break;

  default:
    fprintf(stdout,"readConn: Problem in element type ...\n");
    exit(0);
    break;
  }
}

void
readCoord(FILE *fpi, FILE *fpo)
/***********************************************************************

Routine: readCoord - read in nodal coordinates

Author : Mark A. Christon

Date   : 01/06/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  Real x,y,z;
  int i,nd;
  int Nch,iset;
  char cbuf[MAXCHR];
  char ctok[10];

  /* Parse coordinates */
  fprintf(fpo,"\t Nodal Coordinates\n  ");
  fprintf(fpo,"\t =================\n\n");
  switch (mesh.Ndim) {

  case ONED:
    fprintf(fpo,"\tNode \t x-coord.\n");
    fprintf(fpo,"\t---- \t----------\n");
    for (i=0;i<mesh.Nnp;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%d%g",&nd,&x);
          nd-=1;
          mesh.x[0][nd]=x;
          fprintf(fpo,"\t%4d \t%10.4e\n",nd,mesh.x[0][nd]);
	  iset=1;
	}
#ifdef DEBUG
	else {
	  fprintf(stdout,"%s\n",cbuf);
	}
#endif
      }
    }
    break;

  case TWOD:
    fprintf(fpo,"\tNode \t x-coord. \t y-coord.\n");
    fprintf(fpo,"\t---- \t----------\t---------\n");
    for (i=0;i<mesh.Nnp;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%d%g%g",&nd,&x,&y);
          nd-=1;
          mesh.x[0][nd]=x;
          mesh.x[1][nd]=y;
          fprintf(fpo,"\t%4d \t%10.4e \t%10.4e\n",
		  nd,mesh.x[0][nd],mesh.x[1][nd]);
	  iset=1;
	}
      }
    }
    break;

  case THREED:
    fprintf(fpo,"\tNode \t x-coord. \t y-coord. \t z-coord. \n");
    fprintf(fpo,"\t---- \t----------\t--------- \t----------\n");
    for (i=0;i<mesh.Nnp;i++) {
      iset=0;
      while (iset==0) {
	if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	  sscanf(cbuf,"%d%g%g%g",&nd,&x,&y,&z);
          nd-=1;
          mesh.x[0][nd]=x;
          mesh.x[1][nd]=y;
          mesh.x[2][nd]=z;
          fprintf(fpo,"\t%4d \t%10.4e \t%10.4e \t%10.4e\n",
		  nd,mesh.x[0][nd],mesh.x[1][nd],mesh.x[2][nd]);
	  iset=1;
	}
      }
    }
    break;

  default:
    fprintf(stdout,"ReadCoord: Problem in element type ...\n");
    exit(0);
    break;
  }
}

void
echoHeader(FILE *fp)
/***********************************************************************

Routine: echoHeader - echo the mesh header information

Author : Mark A. Christon

Date   : 01/06/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  fprintf(stdout,"%s\n\n",mesh.title);
  fprintf(stdout,"\tNdim     : %10d\n",mesh.Ndim);
  fprintf(stdout,"\tNnp      : %10d\n",mesh.Nnp );
  fprintf(fp    ,"%s\n\n",mesh.title);
  fprintf(fp    ,"\tNdim     : %10d\n",mesh.Ndim);
  fprintf(fp    ,"\tNnp      : %10d\n",mesh.Nnp );

  switch (mesh.Ndim) {

    case ONED:
      fprintf(stdout,"\tNel      : %10d\n",mesh.Nel2 );
      fprintf(stdout,"\tNmat     : %10d\n",mesh.Nmat2);
      fprintf(stdout,"\tNnpe     : %10d\n",mesh.Nnpe2);
      fprintf(fp    ,"\tNel      : %10d\n",mesh.Nel2 );
      fprintf(fp    ,"\tNmat     : %10d\n",mesh.Nmat2);
      fprintf(fp    ,"\tNnpe     : %10d\n",mesh.Nnpe2);
      break;

    case TWOD:
      fprintf(stdout,"\tNel      : %10d\n",mesh.Nel4 );
      fprintf(stdout,"\tNmat     : %10d\n",mesh.Nmat4);
      fprintf(stdout,"\tNnpe     : %10d\n",mesh.Nnpe4);
      fprintf(fp    ,"\tNel      : %10d\n",mesh.Nel4 );
      fprintf(fp    ,"\tNmat     : %10d\n",mesh.Nmat4);
      fprintf(fp    ,"\tNnpe     : %10d\n",mesh.Nnpe4);
      break;

    case THREED:
      fprintf(stdout,"\tNel      : %10d\n",mesh.Nel8 );
      fprintf(stdout,"\tNmat     : %10d\n",mesh.Nmat8);
      fprintf(stdout,"\tNnpe     : %10d\n",mesh.Nnpe8);
      fprintf(fp    ,"\tNel      : %10d\n",mesh.Nel8 );
      fprintf(fp    ,"\tNmat     : %10d\n",mesh.Nmat8);
      fprintf(fp    ,"\tNnpe     : %10d\n",mesh.Nnpe8);
      break;

    default:
      fprintf(stdout,"\n\tProblem in element type ...\n");
      fprintf(fp    ,"\n\tProblem in element type ...\n");
      exit(0);
      break;
  }

  fprintf(stdout,"\tNnd_sets : %10d\n  ",mesh.Nnd_sets);
  fprintf(stdout,"\tNsd_sets : %10d\n\n",mesh.Nsd_sets);
  fprintf(fp    ,"\tNnd_sets : %10d\n  ",mesh.Nnd_sets);
  fprintf(fp    ,"\tNsd_sets : %10d\n\n",mesh.Nsd_sets);

}

void
ReadMeshNodeset( FILE *fpi, FILE *fpo )
/***********************************************************************

Routine: ReadMeshNodeset - read the nodesets for ascii mesh file

Author : Mark A. Christon

Date   : 04/05/1998

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int  i, j, isread;
  int  ijunk;
  int  Nnd_sets;
  char cbuf[MAXCHR];

  Nnd_sets=mesh.Nnd_sets;

  bzero(ndset,sizeof(Nodeset)*Nnd_sets);

  /* Part 1 of Nodeset data */
  isread=0;
  while (!isread) {
    if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
      sscanf(cbuf,"%d",&ijunk);
      isread=1;
    }
  } 
  if(ijunk!=Nnd_sets) {
    fprintf(stdout,"Error reading nodesets ...\n");
    fprintf(stdout,"Nnd_sets in header: %d\n",Nnd_sets);
    fprintf(stdout,"Nnd_sets in footer: %d\n",ijunk   );
    exit(0);
  }

  /* Part 2 of Nodeset data */
  i=0;
  while (i<Nnd_sets) {
    while((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
      sscanf(cbuf,"%d%d",&ndset[i].id,&ndset[i].Nnp);
      i++;
    }
  }

  /* Part 3 of Nodeset data */
  /* Allocate & read node set lists */
  for (i=0;i<Nnd_sets;i++) {
    ndset[i].node_list=(int *) malloc(ndset[i].Nnp*sizeof(int)); 
    j=0;
    while (j<ndset[i].Nnp) {
      if((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	sscanf(cbuf,"%d%d",&ijunk, &ndset[i].node_list[j] );
	j++;
      }
    }
  }  /* End of allocate & read node set lists */

  /* Now echo the data to fpo */
  fprintf(fpo,"\n");
  fprintf(fpo,"\tNodeset Parameters:\n\t===================\n");
  fprintf(fpo,"\n\tNumber of Node Sets :%d\n",Nnd_sets);
  fprintf(fpo,"\n");
  fprintf(fpo,"\tNode Set#  Node Set Id\n");
  fprintf(fpo,"\t=========  ===========\n");
  for (i=0;i<Nnd_sets;i++) {
    fprintf(fpo,"\t%8d         %4d\n", i+1, ndset[i].id);
  }
  fprintf(fpo,"\n");

  for (i=0;i<Nnd_sets;i++) {
    fprintf(fpo,"\n\tNodeset - %d\n",ndset[i].id);
    fprintf(fpo,"\tNodes\n\t=====\n");
    for (j=0;j<ndset[i].Nnp;j++) {
       fprintf(fpo,"\t%5d\n",ndset[i].node_list[j]);
    } 
  }

}

void
ReadMeshSideset( FILE *fpi, FILE *fpo )
/***********************************************************************

Routine: ReadMeshSideset - read the sidesets for ascii mesh file

Author : Mark A. Christon

Date   : 04/05/1998

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int  i,j,isread;
  int  ijunk;
  int  Nsd_sets;
  char cbuf[MAXCHR];

  Nsd_sets=mesh.Nsd_sets;
  bzero(sdset,sizeof(Sideset)*Nsd_sets);

  /* Part 1 of Sideset data */
  isread=0;
  while (!isread) {
    if((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!= 0)) {
      sscanf( cbuf, "%d", &ijunk );
      isread = 1;
    }
  } 
  if (ijunk!=Nsd_sets) {
    fprintf( stdout, "Error reading sidesets ...\n" );
    fprintf( stdout, "Nsd_sets in header: %d\n", Nsd_sets );
    fprintf( stdout, "Nsd_sets in footer: %d\n", ijunk    );
    exit(0);
  }

  /* Part 2 of Sideset data */
  i = 0;
  while (i<Nsd_sets) {
    while ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
      sscanf(cbuf,"%d%d",&sdset[i].id,&sdset[i].Nel);
      i++;
    }
  }

  /* Part 3 of Sideset data */
  /* Allocate & read side set lists */
  for (i=0;i<Nsd_sets;i++) {
    sdset[i].el  =(int *) malloc(sdset[i].Nel*sizeof(int));
    sdset[i].side=(int *) malloc(sdset[i].Nel*sizeof(int));
    j = 0;
    while (j<sdset[i].Nel) {
      if ((getLine(fpi,cbuf,MAXCHR)>1) && (iscomment(cbuf)!=0)) {
	sscanf(cbuf,"%d%d",&sdset[i].el[j],&sdset[i].side[j]);
	j++;
      }
    }
  }  /* End of allocate & read node set lists */

  /* Now echo the data to fpo */
  fprintf(fpo,"\n");
  fprintf(fpo,"\tSideset Parameters:\n\t================\n");
  fprintf(fpo,"\n\t   Number of Side Sets : %10d\n",Nsd_sets);
  fprintf(fpo,"\n");
  fprintf(fpo,"\t   Side Set#  Side Set Id\n");
  fprintf(fpo,"\t   =========  ===========\n");
  for (i=0;i<mesh.Nsd_sets;i++) {
    fprintf(fpo,"\t    %8d         %4d       %d \n",i+1,sdset[i].id);
  }
  fprintf(fpo, "\n" );
}

int
getLine( FILE *fp, 	/* File pointer for input                    */
         char *s, 	/* Character buffer to be loaded             */
         int  lim  	/* Buffer size in characters                 */)
/***********************************************************************

Routine: getLine - read a line of characters from a file

   Returns an integer count of characters read.

Author : Mark A. Christon

Date   : c 1990

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/

/* This is the getLine function from Kernighan and Ritchie, pg. 67,
 * modified to get characters from the file pointed to by fp.
 */

{
  int c, i;

  i = 0;
  while( --lim > 0 && (c=getc(fp)) != EOF && c != '\n' ) s[i++] = c;

  /* Handle the terminating characters */
  if( c == '\n') s[i++] = '\0';
  
  return( i ); /* Return a count of characters read */

}

int
iscomment( char *cdat )
/***********************************************************************

Routine: iscomment - check an input line for a comment

Author : Mark A. Christon

Date   : 03/1998

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
   int yesno;
   yesno = 1; /* Default is no not a comment */
   if ((strncasecmp(cdat,"#",1 )==0)|| 
       (strncasecmp(cdat,"*",1 )==0)||
       (strncasecmp(cdat,"$",1 )==0)  ) yesno=0;
   return(yesno);
}

int
bandwidth()
/***********************************************************************

Routine: compute the 1/2-bandwidth

Author : Mark A. Christon

Date   : 01/08/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int i, j, ib, ibe;

  ib=0;
  switch (mesh.Ndim) {
    case ONED:
      for (i=0;i<mesh.Nel2;i++) {
        ibe=abs(oned[i].ix[0]-oned[i].ix[1]);
        if (ibe>ib) ib=ibe;
      }
      break;

    case TWOD:
      break;

    case THREED:
      break;

    default:
      fprintf(stdout,"bandwidth: Problem in element type ...\n");
      exit(0);
      break;
  }
  return(ib+1);
}

void
formknf( const Real T,     /* Wire tension */
	 const int ltype,  /* load type =0 constant, =1 linear, =2 split */
	 Real *k,          /* stiffness */
	 Real *f )         /* rhs (force) */
/***********************************************************************

Routine: formknf - form & assemble the K and F matrices (UDU format)

Author : Mark A. Christon

Date   : 02/08/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int i, j, n1, n2, nnp, ibo;
  Real *x;
  Real le, kij, f1, f2, w0, w1, w2;

  nnp=mesh.Nnp;
  x=mesh.x[0];
  /* Loop over elements to assemble */
  for (i=0;i<mesh.Nel2;i++) {
    n1=oned[i].ix[0];
    n2=oned[i].ix[1];
    ibo=abs(n2-n1);
    n1-=1;
    n2-=1;
    /* Setup stiffness first */
    le=fabs(x[n2]-x[n1]);
    kij=T/le;
    /* Only assemble the diagonal & upper-diagonal terms */
    k[n1]+=kij;
    k[n2]+=kij;
    if (n1>n2) {
      k[n2+nnp*ibo]-=kij;
    }
    else {
      k[n1+nnp*ibo]-=kij;
    }

    /* Now compute load */
    switch (ltype) {
      case 0:
        w0=10.0;
	f1=f2=0.5*w0*le;    
	f[n1]+=f1; 
	f[n2]+=f2; 
        break;

      case 1:
        /* linear loading goes here */

        /* BLEEPED CODE GOES HERE  */
	break;

      case 2:
        /* split constant-linear loading goes here */

        /* BLEEPED CODE GOES HERE  */
	break;

      default:
        fprintf(stdout,"formknf: Error in load computation ...\n");
        exit(0);
        break;
    }
  }
}

void
blastbc( Real *k, Real *f ) 
/***********************************************************************

Routine: blastbc - use penalty method to enforce bc's

Author : Mark A. Christon

Date   : 02/08/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int nnp;

  nnp=mesh.Nnp;

  /* all homogeneous EBC's for wire problem */
  k[0]=k[0]*1.0e+10;
  k[nnp-1]=k[nnp-1]*1.0e+10;
  f[0]=0.0;
  f[nnp-1]=0.0;
}

void
write_u( Real *u, FILE *fpo, FILE *fpg)
/***********************************************************************

Routine: write_u - write output and plot files

Author : Mark A. Christon

Date   : 02/08/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int i;

  fprintf(fpo,"\n\tNodal Displacements   ");
  fprintf(fpo,"\n\t===================\n\n");
  fprintf(fpo,"\n\tNode  Displacement\n");
  fprintf(fpo,"  \t====  ============\n");

  fprintf(fpg,"#\n");
  fprintf(fpg,"# Nnp=%4d\n",mesh.Nnp );
  fprintf(fpg,"# Nel=%4d\n",mesh.Nel2);
  fprintf(fpg,"#\n");
  fprintf(fpg,"# Coord.  Displacement\n");
  fprintf(fpg,"#\n");

  for (i=0;i<mesh.Nnp;i++) {
    fprintf(fpo,"\t%4d  %10.4e\n",i,u[i]);
    fprintf(fpg,"%10.4e  %10.4e\n",mesh.x[0][i],u[i]);
  }
}