#include <stdlib.h>

/* zero index version of st; C does not allow for the declaration of array sizes that depend upon
other variables that are passed into a routine.  This cover definition allows one to get around this
difficulty.  stf is a pointer to an array with numnp rows and 6 column which is stored in column
major format.  st(,) then operates just as it would if it were a FORTRAN array */

#define st(a,b) (*(stf+a+(b)*numnp))

/* The definitions for shps and sig are hard wired for 4 Gauss points which they shouldn't be */
void stcn13(int *ix, double shps[4][4][3], double *xsj,
            double sigl[4][6],
            double *dt, double *stf,
            int lint,int nel, int numnp) {

/*
c-------------------------------------------------------------------
c
c      * * F E A P * * A Finite Element Analysis Program
c      Copyright (c) 1984-2000: Robert L. Taylor
c
c      C Language Example Element for FEAP
c      Copyright (c) 2000     : Sanjay Govindjee
c
c      Stress projection routine
c
c      Inputs:
c         ix               - connectivity information
c         shps(lint,nel,3) - shape fncs and derivatives
c         xsj(lint)        - Gauss point jacobians
c         sigl(lint,6)     - Gauss point stresses
c         lint             - number of Gauss points
c         nel              - number of element nodes
c         numnp            - number of nodes in problem
c
c      Output:
c         dt(numnp)        - nodal weights
c         st(numnp,6)      - nodal stresses (accessed throught cover
c                            function and pointer stf)
c
c-------------------------------------------------------------------
*/

      int i, l, ll;

      double  xj;

      /* compute stress projections to nodes */
      for(l = 0 ; l < lint ; l++ ) {

      /* compute projections: int ( sig * shp(i) * darea ) */

        for(i = 0 ; i < nel ; i++ ) {
          ll = abs(ix[i]);
          if(ll != 0) {
            ll -= 1;   /* account for zero based indexing in C by shifting ll back by one location*/
            xj       = xsj[l]*shps[l][i][2];
            dt[ll]   +=  xj;
            st(ll,0) +=  sigl[l][0]*xj;
            st(ll,1) +=  sigl[l][1]*xj;
            st(ll,2) +=  sigl[l][2]*xj;
            st(ll,3) +=  sigl[l][3]*xj;
            st(ll,4) +=  sigl[l][4]*xj;
            st(ll,5) +=  sigl[l][5]*xj;
          }
        }
      }
 

}