/*
 * 
 * ABAQUS FILE OUTPUT reader in C-language
 * 
 * Written by Sami Saarinen, Centre for Scientific Computing (CSC), Finland
 * During Spring'92
 *
 */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include "abaqus.h"

#define MAX(a,b) ( (a) > (b) ? (a) : (b) )
#define MIN(a,b) ( (a) < (b) ? (a) : (b) )

#define PRIVATE static
#define PUBLIC  

typedef union {
  REAL array[513];
  int  jrray[513][PRECISION];
} Array;

PRIVATE Array   data;

#define JRRAY(k)  (data.jrray[k-1][0])
#define ARRAY(k)  (data.array[k-1])

PRIVATE bool    isopen = FALSE;
PRIVATE int     *nw = &JRRAY(1);
#define NW     (*nw)
PRIVATE int     *key= &JRRAY(2);
#define KEY    (*key)
#define NATTR  ( NW - 2 )

#define NODEID  (*nodeid)
#define ELEMID  (*elemid)

PRIVATE int     numnodes  = 0;
PRIVATE int     numelems  = 0;
PRIVATE int     maxnodeid = 0;
PRIVATE int     maxelemid = 0;

/*
 * 
 * Routines with internal linkage i.e. static/PRIVATE to this file 
 *
 */


PRIVATE int   CD_Original(void);
PRIVATE char *Parse_Path(char *file);
PRIVATE int   Open(char *file);
PRIVATE int   Close(void);
PRIVATE int   Rewind(void);
PRIVATE int   Read_to_Array(void);
PRIVATE int   Locate_to_Key(int keynum);
PRIVATE char *Get_Elem_Characteristics(REAL elemtype,
				      Descr *geom,
				      Descr *order,
				      Descr *str_flag,
				      int numnodes_per_elem);
PRIVATE char *Get_String(REAL hollerith, bool strip_flag);
PRIVATE int   Read_Pure_Nodal_Data(NodalData *nd, 
				   int keynum, 
				   bool rewind_flag,
				   bool alloc_flag,
				   int ndatacomp);
PRIVATE void  Zero_Nodal_Data(NodalData *nd);
PRIVATE unsigned int Hash(char *s);
/* PRIVATE unsigned int Lookup(char *s); inlined ...! */
#define HASHSIZE    1013
PRIVATE unsigned int Hashtable[HASHSIZE];


PRIVATE unsigned int Hash(char *s)
/*
 *
 * Calculate hash value (taken directly from K & R) 
 * Enlarge HASHSIZE to larger prime number if gives
 * same hash value with definitely different string s
 *
 */
{
  unsigned int hashval;
  for (hashval=0; *s ; s++) hashval = (*s) + 31 * hashval;
  return hashval % HASHSIZE;
}

/*
PRIVATE unsigned int Lookup(char *s)
{
  return ( Hashtable[Hash(s)] );
}
*/

#define Lookup(s) (Hashtable[Hash((s))])

PRIVATE char *Get_String(REAL hollerith, bool strip_flag)
/*
 *
 * A quick and dirty way to extract string stored in a 64-bit float.
 * Substitutes ASCII NULL to first occurence of SPACE-char if
 * strip_flag is TRUE.
 *
 */
{
  typedef union {
    char s[9];
    REAL h;
  } Share;

  static Share it;
  char *p;

  it.h = hollerith;
  if (strip_flag) {
    it.s[8] = ' ';
    p = it.s;
    while ( *p++ != ' ') ;
    *--p = '\0';
  }
  it.s[8] = '\0';

  return (it.s);
}

PRIVATE int CD_Original()
/*
 * Change Directory to original
 */
{
  static bool first_call = TRUE;
  static char original[BUFSIZ] = "./"; /* BUFSIZ as in stdio.h */
  int retcode = OK;

  if (first_call) {
   getcwd(original,BUFSIZ); 
   first_call = FALSE;
  }
  else {
    if (chdir(original) != 0) {
      perror(original);
      retcode = ERROR;
    }
  }
  return (retcode);
}

PRIVATE char *Parse_Path(char *file)
/*
 * Parses pathname from file. If pathname is not found
 * a "./" substituted (i.e. current directory).
 * Also, pathname is stripped away from *file.
 */
{
  register char *p = NULL;
  static char *pathname = NULL;

  if (pathname) free(pathname);

  if ((p = strrchr(file,'/')) == NULL) {
    pathname = (char *)calloc(3,sizeof(char));
    strcpy(pathname,"./");
  }
  else {
    register char *strippedfile = file;
    register int len=0;
    register char *strippedpath = NULL;

    while ( strippedfile != p ) { /* Yeah, I meant that, a pointer comparison ! */
      len++;
      strippedfile++;
    }
    pathname = (char *)calloc(len+1,sizeof(char));
    strippedfile = file;
    strippedpath = pathname;
    while ( strippedfile != p ) { /* Yeah, I meant that, a pointer comparison ! */
      *strippedpath++ = *strippedfile++;
    }
    *strippedpath = '\0';

    strippedfile = file;
    do {
      *strippedfile++ = *++p;
    } while ( *p );
  }

  return pathname;
}

PRIVATE int   Open(char *file)
/*
 *
 * Open ABAQUS File Output
 *
 */
{
  register int i,len;
  int nru=1, loutf=0, junit=8;
  int lrunit[2];
  FILE *fp;
  char *pathname = NULL;
  static char *filename = NULL;
#ifdef __convex__
  static char *FOR008 = "for008.fil";
  static bool dup2_trick = FALSE;
#endif

  extern struct {
    char name[26];
  } JOBPAR;
  

  if (isopen) {
    /* fprintf(stderr,"File Output is already open\n"); */
    return(OK);
  }
  else
    isopen = FALSE;


  if (!file && !filename) {
    fprintf(stderr,"File Output name wasn't given\n");
    return(ERROR);
  }

  if (file) {
    pathname = Parse_Path(file);

    CD_Original();
    if (chdir(pathname) != 0) {
      perror(file);
      return(ERROR);
    }

    if (filename) free(filename);

    /* 
     *  Enables opening of several FILE OUTPUT files
     *  in a single execution after previous ones have 
     *  been closed first.
     */

    filename = (char *)calloc(strlen(file) + 1, sizeof(char));
    strcpy(filename,file);
#ifdef DEBUG
    fprintf(stderr,">>New filename='%s'\n",filename);
#endif

#ifdef __convex__
    {
      if ( strcmp(filename,"$stdin") == 0 && !dup2_trick) {
	free(filename);
	filename = (char *)calloc(strlen(FOR008) + 1, sizeof(char));
	strcpy(filename,FOR008);
	fp = fopen(filename,"wb");
	if (!fp) {
	  fprintf(stderr,"Cannot create temporary File Output '%s' \n",filename);
	  return(ERROR);
	}
	fclose(fp);
	dup2_trick = TRUE;
      }
    }
#endif

    if ( !(fp = fopen(filename,"rb") ) ) {
      fprintf(stderr,"File Output '%s' not found\n",filename);
      return(ERROR);
    }
    else
      fclose(fp);

    len = sizeof(JOBPAR.name);
    for (i=0 ; i<len ; i++) JOBPAR.name[i] = ' ';
    
    len = strlen(filename);
    len = MIN(len,sizeof(JOBPAR.name) - 1);
    
    for (i=len-1; i>0 ; i--) {
      if (filename[i] == '.') {
	len = i;
	break;
      }
    }
    for (i=0; i<len; i++)  JOBPAR.name[i] = filename[i];
    len = sizeof(JOBPAR.name) - 1;
    JOBPAR.name[len] = '\0';
    
  }

  lrunit[0] = 8;  
  lrunit[1] = 2;
    
  INITPF(&nru,lrunit,&loutf);
  DBRNU(&junit);

#ifdef __convex__
  if (dup2_trick)
  {
    extern int for$getfp_(int *fortran_unit);
    int fp_ptr = for$getfp_(&junit); /* Get FILE pointer of Fortran UNIT #8 */
    FILE *fp = (FILE *)fp_ptr;
    if (!fp) {
      fprintf(stderr,"Error in reading File Output from standard input\n");
      return(ERROR);
    }
    /* Close temporary file for008.fil and initiate reading from stdin */
    dup2(fileno(stdin),fileno(fp)); 
    fp = stdin;
  }
#endif

  isopen = TRUE;
  return(OK);

}

PRIVATE int Close()
/*
 *
 * Marks File Output file to be closed
 *
 */
{
  isopen = FALSE;
  return(OK);
}

PRIVATE int Rewind()
/*
 *
 * Go to beginning of File Output
 *
 */
{
  int rc=0;
/*
 * Rewind thru CALL DBFILE(0,xxx,retcode) seemed not to work neither on Convex
 * or Iris, so it was replaced by 'safer' Close + Open
 *
  static int two=2;
  REAL dummy;
  DBFILE(&two, &dummy, &rc);
 */
  Close();
  Open(NULL);
  return (rc);
}

PRIVATE int Read_to_Array()
/*
 *
 * Read one logical/physical record from File Output
 *
 */
{
  static int zero=0;
  int rc;
  DBFILE(&zero, &ARRAY(1), &rc);
  return (rc);
}

PRIVATE int Locate_to_Key(int keynum)
/*
 *
 * Locate File Output into the beginning of the next occurence 
 * of key number KEY which equals to keynum
 *
 */
{
  int rc;
  do {
    rc = Read_to_Array();
  } while ( rc == 0 && KEY != keynum );
  
  return ( ( rc==0 ) ? keynum : ERROR );
}


PRIVATE char *Get_Elem_Characteristics(REAL elemtype,
				      Descr *geom,
				      Descr *order,
				      Descr *str_flag,
				      int numnodes_per_elem)
/*
 *
 * Contains all ABAQUS-elements currently supported by the interface 
 * First, hash table is set up, containing indices to SE-table.
 * Then each elemtype, converted to string, is in turn hashed to get
 * their index via Lookup. If index > 0, then the data if found at SE[index].
 * If index == 0, then the unknown element characteristics are return.
 *
 */
{
  typedef struct {
    char  *name;
    Descr geom;
    Descr order;
    Descr str_flag;
  } Supported_Elems;

  static Supported_Elems SE[] = {
#include "eltypes.h"
  };

  static bool first_time=TRUE;

  Supported_Elems *sep;
  unsigned int index;

  if (first_time) {
    register int i;
    for (i=0; i<HASHSIZE; i++) 
      Hashtable[i] = 0;
    for (i=0, sep=SE; sep->name ; sep++, i++)
      Hashtable[Hash(sep->name)] = i;
#ifdef DEBUG
    for (i=0; i<HASHSIZE; i++) {
      if ( (index=Hashtable[i]) > 0) {
	sep = &SE[index];
	fprintf(stderr,"Index=%d, Name='%s', Hahsval=%u\n",
		index,sep->name,i);
      }
    }
#endif    
    first_time = FALSE;
  }

  if ( ( index = Lookup(Get_String(elemtype,TRUE)) ) > 0 ) 
    sep = &SE[index];
  else 
    sep = SE;

  *geom = sep->geom;
  *order = sep->order;
  *str_flag = sep->str_flag;
  return (sep->name);
}



PRIVATE void Zero_Nodal_Data(NodalData *nd)
/* 
 * Fast zeroing of data area. Assumes contiguous memory alloc. 
 */
{
  memset(nd->data[0], 0 , nd->ndatacomp * nd->numnodes * sizeof(REAL));
}


PRIVATE int Read_Pure_Nodal_Data(NodalData *nd,
				 int keynum, 
				 bool rewind_flag,
				 bool alloc_flag,
				 int ndatacomp)
/*
 *
 * A slave routine to read most of the nodal data records available on ABAQUS.
 * Reads at most MIN(ndatacomp,NATTR-1) data components per particular node.
 * Assumes that required arrays are already allocated if alloc_flag is FALSE.
 *
 */
{
  register int rc;

  if (!nd) return(ERROR);

  if (!isopen) {
    fprintf(stderr,"ABAQUS File Output wasn't opened\n");
    return(ERROR);
  }

  if (rewind_flag) Rewind();

  if (alloc_flag) rc = ABAQUS_Alloc_NodalData(nd,ndatacomp);
  
  if (rc == ERROR) {
    fprintf(stderr,"Nodal data allocation failed\n");
    return(ERROR);
  }

  Zero_Nodal_Data(nd);

  if (Locate_to_Key(keynum) == keynum) {
    register int j=0;
    register REAL *ndp;
    register int i, nmax;
    int *nodeid = &JRRAY(3);
    do {
      ndp = nd->data[j];
      nd->nid[j] = NODEID;
      nmax = MIN(ndatacomp,NATTR-1);
      for (i=0; i<nmax; i++) ndp[i] = ARRAY(i+4);
      if (nmax < ndatacomp) 
	for (i=nmax; i<ndatacomp; i++) ndp[i] = 0.0;

      rc = Read_to_Array();
      j++;
    } while (rc == 0 && KEY == keynum);
    return(OK);
  }
  else {
    fprintf(stderr,"Nodal data set #%d not found\n",keynum);
    return(ERROR);
  }
}


/*
 * 
 * Routines with external linkage 
 *
 */

PUBLIC int ABAQUS_Debug(FILE *fp, 
			bool rewind_flag_start,
			bool rewind_flag_end
			)
/*
 *
 * Reads File Output from the current or rewound (at the beginning)
 * position and prints out record headers record by record until
 * EOF is reached. 
 * For debugging purposes only.
 *
 */
{
  int nrecords=0;
  int prev_key=-1;

  if (!isopen) {
    fprintf(stderr,"ABAQUS File Output wasn't opened\n");
    return(ERROR);
  }

  if (rewind_flag_start) Rewind();

  fprintf(fp,"RECORD#   Key  NW  Nattr\n\n");
  
  while ( Read_to_Array() == 0 ) {
    nrecords++;
    if (prev_key == KEY)
      fprintf(fp,"%7d%10d%7d\n"  ,nrecords,    NW,NATTR);
    else {
      fprintf(fp,"%7d%6d%4d%7d\n",nrecords,KEY,NW,NATTR);
      prev_key = KEY;
    }
  }

  if (rewind_flag_end) Rewind();

  return (nrecords);
}


PUBLIC int ABAQUS_Open(char *file)
/*
 *
 * Opens File Output file pointed by filename.
 * ABAQUS-library routines assume (under Unix) that Fortran
 * COMMON/JOBPAR/ storing string of 25 characters, contains
 * filename W/O extension '.fil'. That is stripped away in the
 * beginning (after the file existence is tested).
 *
 */
{

  if (Open(file) > 0) {
  
    Rewind();
    if (Locate_to_Key(2001) != 2001) {
      fprintf(stderr,"Warning: Element & Nodal data not written yet\n");
      Rewind();
      return(NODATA);
    }

    Rewind();
    if (Locate_to_Key(1921) == 1921) {
      register int i;
      register REAL *version = &ARRAY(3);
      register int *nelems = &JRRAY(7);
      register int *nnodes = &JRRAY(8);
      
      fprintf(stderr,"Ran by ABAQUS version %8s\nAnalysis Date: ",
	      Get_String(*version++,TRUE));
      for (i=1; i<4; i++, version++)
	fprintf(stderr,"%8s",Get_String(*version,FALSE));
      fprintf(stderr,"\n");
      fprintf(stderr,"# of    nodes in model: %6d\n",*nnodes);
      fprintf(stderr,"# of elements in model: %6d\n",*nelems);
      
      numnodes = *nnodes;
      numelems = *nelems;
      
      Rewind();
      if (Locate_to_Key(1922) == 1922) {
	register REAL *header = &ARRAY(3);
	
	for (i=0; i<10; i++, header++)
	  fprintf(stderr,"%8s",Get_String(*header,FALSE));
	fprintf(stderr,"\n");
      }
      else
	fprintf(stderr,"Warning: Can't locate ABAQUS header record\n");
      
      Rewind();
      return(OK);
    }
    else {
      fprintf(stderr,"Can't locate ABAQUS version & model size data record\n");
      fprintf(stderr,"File Output possibly closed or corrupted\n");
      isopen = FALSE;
      return(ERROR);
    }
  }
  else {
    fprintf(stderr,"Open to file '%s' failed\n",(file)? file : "");
    return(ERROR);
  }

}


int ABAQUS_Rewind()
/*
 *
 * Rewinds ABAQUS File Output
 *
 */
{
  return (Rewind());
}

int ABAQUS_Close()
/*
 *
 * Marks File Output file to be closed
 *
 */
{
  return(Close());
}



PUBLIC Node *ABAQUS_Read_Nodes()
/*
 *
 * Reads nodal coordinate data
 * Assumes all specific data in subsequent records
 *
 */
{
  Node *p = NULL; /* Pointer to Node-structure array to be allocated */

  if (!isopen) {
    fprintf(stderr,"ABAQUS File Output wasn't opened\n");
    return(NULL);
  }

  Rewind();

  /* Allocate room for nodal coordinates */

  if ( Locate_to_Key(1901) == 1901 ) {
    p = (Node *) calloc(numnodes, sizeof(Node));
    if (!p) {
      fprintf(stderr,"Can't allocate memory for nodal coordinates\n");
      fprintf(stderr,"Storage needed was %d bytes\n",
	      numnodes * sizeof(Node));
      return(NULL);
    }
  }
  else {
    fprintf(stderr,"Can't find nodal coordinate records\n");
    return(NULL);
  }


  /* Read actual data */

  {
    register int i=0;
    register int *nodeid = &JRRAY(3);
    register int rc;
    do {
      maxnodeid = MAX(maxnodeid,NODEID);
      p[i].nid = NODEID;
      p[i].x = ARRAY(4);
      p[i].y = ARRAY(5);
      p[i].z = ARRAY(6);
      rc = Read_to_Array();
      i++;
    } while ( rc == 0 && KEY == 1901 );
  }

  return (p);
}



PUBLIC Element *ABAQUS_Read_Elements()
/*
 *
 * Reads element topology data in 2 passes:
 * Pass #1: Determine total length of nodelist arrays summed over all elements
 * Pass #2: Allocate nodelist arrays and fill in data
 * Assumes all specific data in subsequent records
 *
 */
{
  Element *p = NULL;     /* Pointer to element data structure array to be allocated */
  int *enodelist = NULL; /* Pointer to full nodal connectivity list                 */

  if (!isopen) {
    fprintf(stderr,"ABAQUS File Output wasn't opened\n");
    return(NULL);
  }


  /* Define storage needed to save nodal connectivity list  */

  Rewind();

  if ( Locate_to_Key(1900) == 1900 ) {
    register int storage=0;
    register int rc;
    do {  
      storage += NW - 4;
      rc = Read_to_Array();
    } while ( rc == 0 && KEY == 1900 );


    p = (Element *)calloc(numelems, sizeof(Element));
    if (!p) {
      fprintf(stderr,"Can't allocate memory for element topology array\n");
      fprintf(stderr,"Storage needed was %d bytes\n",
	      numelems * sizeof(Element));
      return(NULL);
    }

    enodelist = (int *) calloc(storage, sizeof(int));
    if (!enodelist) {
      fprintf(stderr,"Can't allocate memory for full nodal connectivity list\n");
      fprintf(stderr,"Storage needed was %d bytes\n",
	      storage * sizeof(int));
      return(NULL);
    }

  }
  else {
    fprintf(stderr,"Can't find element topology records\n");
    return(NULL);
  }


  /* Read actual data */

  Rewind();
  Locate_to_Key(1900);

  {
    register int i=0, j;
    register int *elemid = &JRRAY(3);
    register REAL *elemtype = &ARRAY(4);
    register int rc;

    do {
      maxelemid = MAX(ELEMID,maxelemid);
      p[i].eid = ELEMID;
      p[i].numnodes = NW - 4;
      p[i].name = Get_Elem_Characteristics(*elemtype,
					   &p[i].geom,
					   &p[i].order,
					   &p[i].str_flag,
					   p[i].numnodes);
      p[i].nodelist = enodelist;
      for ( j = 0; j < p[i].numnodes; j++) enodelist[j] = JRRAY(j+5);
      enodelist += p[i].numnodes;
      rc = Read_to_Array();
      i++;
    } while ( rc == 0 && KEY == 1900 );
    
  }
  return (p);
}


PUBLIC int ABAQUS_Alloc_NodalData(NodalData *nd, 
				  int ndatacomp)
/*
 *
 * Allocates room for nodal data set, ndatacomp values per node
 * This MUST BE called before use of particular nodal data set begins
 *
 */
{
  register int storage = numnodes * ndatacomp;
  register REAL *ndp;
  
  if (!nd) return(ERROR);

  nd->data = (REAL **) calloc(numnodes, sizeof(REAL *));
  if (!nd->data) {
    fprintf(stderr,"Can't allocate storage for data set data pointers\n");
    fprintf(stderr,"Storage needed was %d bytes\n",
	    numnodes * sizeof(REAL *));
    return(ERROR);
  }

  nd->nid = (int *) calloc(numnodes,sizeof(int));
  if (!nd->nid) {
    fprintf(stderr,"Can't allocate storage for data set node id array\n");
    fprintf(stderr,"Storage needed was %d bytes\n",
	    numnodes * sizeof(int));
    return(ERROR);
  }
    
  ndp = (REAL *) calloc(storage, sizeof(REAL)); /* Contiguous memory alloc */

  if (!ndp) {
    fprintf(stderr,"Can't allocate memory for nodal data set data arrays\n");
    fprintf(stderr,"Storage needed was %d bytes\n",
	    storage * sizeof(REAL));
    return(ERROR);
  }

  /* Setup pointers correctly */

  {
    register int i;
    for ( i = 0; i < numnodes; i++, ndp += ndatacomp) 
      nd->data[i] = ndp;
  }

  nd->ndatacomp = ndatacomp;
  nd->numnodes = numnodes;

  return (OK);
}

PUBLIC int ABAQUS_Free_NodalData(NodalData *nd)
/*
 *
 * Free NodalData storage allocated previously by ABAQUS_Alloc_NodalData
 *
 */
{
  if (!nd) return(ERROR);

  free(nd->data[0]);
  free(nd->data);

  return(OK);
}


PUBLIC int ABAQUS_Read_Disp(NodalData *nd)
/*
 *
 * Reads displacement data set (key #101)
 *
 */
{
  return ( Read_Pure_Nodal_Data(nd,101,
				FALSE,
				FALSE,
				nd->ndatacomp) );
}


PUBLIC int ABAQUS_Read_Temp(NodalData *nd)
/*
 *
 * Reads temperature data set (key #201)
 *
 */
{
  return ( Read_Pure_Nodal_Data(nd,201,
				FALSE,
				FALSE,
				nd->ndatacomp) );
}


PUBLIC int ABAQUS_Locate_Time_Step(NodalData *nd, int step, int inc)
/*
 *
 * Locates FILE OUTPUT to a specific time STEP at its increment INC
 *
 */
{
  register int retcode;

  do {
    if ( (retcode = ABAQUS_Next_Time_Step(nd)) > 0 ) {
      if ( nd->step == step && nd->inc == inc ) break;
    }
  } while ( retcode != ERROR );

  if (retcode == ERROR) {
    Rewind();
    do {
      if ( (retcode = ABAQUS_Next_Time_Step(nd)) > 0 ) {
	if ( nd->step == step && nd->inc == inc ) break;
      }
    } while ( retcode != ERROR );
  }
  return (retcode);
}

PUBLIC int ABAQUS_Next_Time_Step(NodalData *nd)
/*
 *
 * Locates next time step start in File Output
 * but DOESN'T check whether step has been written
 * to FILE OUTPUT completely.
 *
 */
{
  if (!nd) return(ERROR);

  if (Locate_to_Key(2000) == 2000) {
    REAL *time = &ARRAY(3);
    int *proc = &JRRAY(7);
    int *step = &JRRAY(8);
    int *inc = &JRRAY(9);
#ifdef DEBUG
    fprintf(stderr,
	    "Procedure #%d, Time=%.5f, Step=%d, Inc=%d\n", 
	    *proc, *time, *step, *inc);
#endif
    nd->step = *step;
    nd->inc = *inc;
    nd->time = *time;
    return (OK);
  }
  else
    return(ERROR);
}



PUBLIC int ABAQUS_Time_Data(STime *timeinc)
/*
 *
 * Reads all currently available time increment data to timeinc in 2 passes:
 * Pass #1: Determine number of incs of FULLY COMPLETED steps
 * Pass #2: Allocate and fill arrays
 *
 */
{
  REAL *time = &ARRAY(3);
  int *proc = &JRRAY(7);
  int *step = &JRRAY(8);
  int *inc = &JRRAY(9);

  if (!timeinc) return(ERROR);

  Rewind();

  timeinc->numincs = 0;

  while ( Locate_to_Key(2000) == 2000 ) { /* Pass #1 */
    if ( Locate_to_Key(2001) == 2001 ) timeinc->numincs++;
    else break;
  }

  if (timeinc->numincs == 0) return(ERROR);

  timeinc->step = (int *)calloc(timeinc->numincs,sizeof(int));
  timeinc->inc = (int *)calloc(timeinc->numincs,sizeof(int));
  timeinc->time = (REAL *)calloc(timeinc->numincs,sizeof(REAL));

  Rewind();
  {
    register int i=0;
    while ( Locate_to_Key(2000) == 2000 ) { /* Pass #2 */
      fprintf(stderr,"Fully completed Step #%d, Inc #%d, ",*step,*inc);
      fprintf(stderr,"Time=%.5f\n", *time);
      timeinc->step[i] = *step;
      timeinc->inc[i] = *inc;
      timeinc->time[i]= *time;
      if ( ++i >= timeinc->numincs) break;
    }
  }
  return(OK);
}



PUBLIC int ABAQUS_Read_Model(Model *model)
/*
 *
 * Reads node and element data using slave routines
 *
 */
{
  if (!model) return(ERROR);
  if (!(model->node = ABAQUS_Read_Nodes()))    return(ERROR);
  if (!(model->elem = ABAQUS_Read_Elements())) return(ERROR);
  model->maxnodeid = maxnodeid;
  model->maxelemid = maxelemid;
  model->numnodes = numnodes;
  model->numelems = numelems;
  return(OK);
}



int ABAQUS_Read_Stress(NodalData *str,
		       int str_flag, int str_type, int str_layer)
/*
 *
 * Reads ABAQUS Stress Data Averaged at Nodes
 *
 * str:      ptr to stress result data structure allocated before.
 * str_flag: stresses to be read; one of 
 *           {STR_SOLID, STR_SHELL, STR_BEAM} meaning
 *           stresses on solid, shell and beam elements respectively
 *           Only  STR_SOLID &  STR_SHELL are implemented.
 * str_type: Stress type to be processed; one of
 *           {STR_NORMAL, STR_VON_MISES, STR_MAX_SHEAR} meaning
 *           regular stress-tensor, Von Mises- and Maximum Shear equivalent
 *           stresses.
 * str_layer:If str_flag == STR_SHELL, then results for particular shell
 *           layer is searched. If set to STR_BOTTOM, then shell surface #1
 *           is assumed, if set to STR_TOP (a big number), the top surface
 *           for each shell is assumed.
 *
 */
{
  int keynum;

  if (!str) return(ERROR);

#ifdef DEBUG
  fprintf(stderr,"ABAQUS_Read_Stress\n");
  fprintf(stderr,"ndc=%d flag=%d type=%d layer=%d\n",
	  str->ndatacomp,str_flag,str_type,str_layer);
#endif

  if (str->ndatacomp < 1) return(ERROR); 


  Zero_Nodal_Data(str);

  switch (str_type) {

  case STR_VON_MISES:
  case STR_MAX_SHEAR:
    keynum = 12;
    break;
  case STR_NORMAL:
    keynum = 11;
    break;
  default:
    return(ERROR);

  }

  if (str_flag != STR_SOLID && str_flag != STR_SHELL) return(ERROR);

  if (str_flag == STR_SOLID && str_layer > 0) return(ERROR);

  if (Locate_to_Key(1) == 1) {
  
    register int j = -1;
    register int prev_nodeid =-1;
    register int rc;
    int attr1;
    int *nodeid = &attr1;
    int attr2;
    int layer;
    int attr4;


    do {
    /*
     * If record key == 1 contains nodal averaged values, then
     * it MUST contain following data in its attributes:
     *
     * Attr#   Data description:
     *
     *   1     Node id
     *   2     0 (= zero)
     *   3     0 if STR_SOLID, layer number if STR_SHELL
     *   4     4
     *   5     Ignore!
     *   6     Number of direct stresses at a point (NDI)
     *   7     Number of shear stresses at a point (NSHR)
     *   8     Ignore!
     *   9     Ignore!
     */
  
      attr1 = JRRAY(3);
      attr2 = JRRAY(4);
      layer  = JRRAY(5);
      attr4 = JRRAY(6);
      
      if (attr2 != 0 || attr4 != 4) goto Next;

      if (str_layer > 0 && str_flag == STR_SHELL) {
	if (str_layer != STR_TOP) {
	  if (layer != str_layer) goto Next;
	}
      }

      if (layer != 0 && str_flag == STR_SOLID) goto Next;

      if (Locate_to_Key(keynum) == keynum) {
	register REAL *ndp;
	register int i, nmax;

	if (prev_nodeid != NODEID) {
	  prev_nodeid = NODEID;
	  j++;
	}
	
	ndp = str->data[j];
	str->nid[j] = NODEID;
	nmax = MIN(str->ndatacomp,NATTR);
	
	switch (str_type) {

	case STR_NORMAL:
	  for (i=0; i<nmax; i++) ndp[i] = ARRAY(i+3);
	  if (nmax < str->ndatacomp) 
	    for (i=nmax; i<str->ndatacomp; i++) ndp[i] = 0.0;
	  break;
	case STR_VON_MISES:
	  ndp[0] = ARRAY(3); 
	  break;
	case STR_MAX_SHEAR:
	  ndp[0] = ARRAY(4); 
	  break;
	default:
	  break;

	} /* switch */
	
      } /* if */

    Next:
      
      do {
	rc = Read_to_Array();
      } while (rc == 0 && (KEY != 1 && KEY != 2001));
      
    } while (rc == 0 && KEY == 1 && j < str->numnodes-1); /* do-while */

  } /* if */

  else {
    fprintf(stderr,"Stress data set #%d not found\n",keynum);
    return(ERROR);
  }


  return(OK);
}

