/* $Header: /usr/home/sanjay/Feap/Elem/Ec/RCS/stif13.c,v 1.3 2000/01/28 20:55:02 sanjay Exp $ */

#include <stdio.h>

/* define s(nst,nst) as a zero index array (column major)*/
#define s(a,b) (*(sf+(a)+((b)*nst)))
#define xl(a,b) (*(xlf+(a)+((b)*ndm)))

struct {
    double hr[1];
    int mr[1000];
} _BLNK__;

struct {
  int numnp,numel,nummat,nen,neq,ipr;
} cdata_;

struct {
    char o[4], head[80];
} bdata_;


struct {
  double dm;
  int    n,ma,mct,iel,nel;
} eldata_;

struct {
  double tt[100];
} elplot_;

struct {
  int    nh1, nh2, nh3;
} hdata_;

struct {
    int ior, iow;
} iofile_;

struct {
    int    nph, ner;
    double erav;
} prstrs_;


struct {
   double ttim,dt,c1,c2,c3,c4,c5;
} tdata_;

static char stress_header[] = "\nElement  Material  Stress-11  Stress-22\
  Stress-33  Stress-12  Stress-23  Stress-13\nx-coord  y-coord\
   Strain-11  Strain-22  Strain-33  Strain-12  Strain-23  Strain-13\n";

static char stress_format[] ="%7d%10d%11.3e%11.3e%11.3e%11.3e%11.3e%11.3e\n\
%7.4f%9.4f%12.3e%11.3e%11.3e%11.3e%11.3e%11.3e\n";

void stif13(double (* d), double (* ul), double (* xlf),
             int    (* ix),double (* tl), double (* sf ),
             double (* p ), 
             int ndf, int ndm,
             int nst, int isw) {

/*
c-------------------------------------------------------
c      * * F E A P * * A Finite Element Analysis Program
c      Copyright (c) 1984-1994: Robert L. Taylor
c
c      C Language Example Element for FEAP
c      Copyright (c) 2000     : Sanjay Govindjee
c
c      Input:
c             d, ul, xlf, ix, tl, ndf, ndm, nst, isw
c
c      Output: (various)
c             p, s
c
c-------------------------------------------------------
*/

	double sg[4][3];      /* swapped order from fortran */
	double shps[4][4][3]; /* swapped order from fortran */
        double ds[6][6], sigav[6], sigl[4][6], xsj[4];
        double sig[6], eps[4][6], deps[4][6], dvol;
        double bbd[4][2], epsav[6], xx1, xx2;

        int  tmp1,tmp2,tmp3,FALSE=0;
        int  nhv,lint,l,hpl;

        int  j,j1,i,i1,jj;

        char iow_str[1600];
        int  len;


        switch(isw) {

          case 3:
          case 6:
          case 4:
          case 8:


/*      
c          Integration points and weights 
*/
           tmp1 = 2; /* Gauss points per direction */
           lint = 4; /* Gauss points in element */
           int2d_(&tmp1,&lint,*sg);


/*      
c          Get shape functions and jacobians
c          And linear strains
c          Note the pointer-pointer arithmetic on sg and shps
c          so be careful if you need to modify these calls
*/
           for(l=0; l<lint; l++) 
             shp2d_(*(sg+l), xlf, **(shps+l),    
                   xsj+l,&ndm,&eldata_.nel,ix,&FALSE);

           
           for(l=0 ; l < lint ; l++) 
             kine13(**(shps+l),ul,ndf,eldata_.nel,*(eps+l),*(deps+l)); 

           nhv = d[10];

           switch(isw) {
             case 3:
             case 6:

               hpl = 0;

               for(l=0; l<lint; l++) {

                 modl13(d,*(eps+l), *(deps+l), tdata_.dt, 
                          &(_BLNK__.hr[hdata_.nh1+hpl-1]), 
                          &(_BLNK__.hr[hdata_.nh2+hpl-1]), sig, *ds);


                 /* Store stresses for tplot */
                 j1 = 6*l;
                 for(j=0 ; j<6 ; j++) elplot_.tt[j + j1] = sig[j];


                 /* Multiply stress and tangent by the volume element */
                 dvol = xsj[l]*sg[l][2];
                 for(i=0 ; i<6 ; i++) {
                  sig[i] = sig[i]*dvol;
                  for(j=0 ; j<6 ; j++) ds[i][j] = ds[i][j]*dvol;
                 }
                

                 i1 = 0;
                 for(i=0 ; i<eldata_.nel ; i++) {
                  p[i1]   = p[i1]   - shps[l][i][0]*sig[0]
                                    - shps[l][i][1]*sig[3];
                  p[i1+1] = p[i1+1] - shps[l][i][0]*sig[3]
                                    - shps[l][i][1]*sig[1];
                  i1 = i1 + ndf;
                 }


                 if(isw==3) {
                  i1 = 0;
                  for(i=0 ; i<eldata_.nel ; i++) {

                   /* Compute bmat-t * ds * dvol  */

                   for(jj = 0; jj < 4; jj++) {

                    bbd[jj][0] = shps[l][i][0]*ds[0][jj]
                               + shps[l][i][1]*ds[3][jj];

                    bbd[jj][1] = shps[l][i][0]*ds[3][jj]
                               + shps[l][i][1]*ds[1][jj];

                   }

                   /* Compute tangent stiffness */

                   j1 = 0;
                   for(j = 0;j< eldata_.nel; j++) {

                    for(jj=0;jj<2;jj++) {

                     s(i1+jj,j1  ) = s(i1+jj,j1  ) + bbd[0][jj]*shps[l][j][0]
                                                   + bbd[3][jj]*shps[l][j][1];

                     s(i1+jj,j1+1) = s(i1+jj,j1+1) + bbd[3][jj]*shps[l][j][0]
                                                   + bbd[1][jj]*shps[l][j][1];
                    }

                    j1 = j1 + ndf;
                   } 

                   i1 = i1 + ndf;
                  }

                }

                hpl = hpl + nhv;
               }
               break;

             case 4:
             case 8:
               hpl = 0;
               for(i=0; i<6 ; i++) sigav[i] = 0.0;
               for(i=0; i<6 ; i++) epsav[i] = 0.0;
               for(l=0; l<lint; l++) {
                 modl13(d,*(eps+l), *(deps+l), tdata_.dt, 
                          &(_BLNK__.hr[hdata_.nh1+hpl-1]), 
                          &(_BLNK__.hr[hdata_.nh2+hpl-1]), *(sigl+l), *ds);

                 for(i=0; i<6 ; i++) {
                   sigav[i] += sigl[l][i]*0.25;
                   epsav[i] += eps[l][i]*0.25;
                 }
                 hpl = hpl + nhv;
               }
               if(isw==4) { 

                 xx1 = 0.0;
                 xx2 = 0.0;
                 for(i=0;i<eldata_.nel;i++) {
                     xx1 += xl(1,i)*0.25;
                     xx2 += xl(2,i)*0.25;
                 }

                 eldata_.mct =  eldata_.mct - 2;
                 if(eldata_.mct <= 0) {

                   tmp1 = 4; fortiow_(bdata_.o,&tmp1);
                   tmp1 = 80; fortiow_(bdata_.head,&tmp1);

                   len = sprintf(iow_str,stress_header);
                   fortiow_(iow_str,&len);

                   if(iofile_.ior < 0) {

                     for(i=0;i<4;i++) printf("%c",bdata_.o[i]); 
                     for(i=0;i<80;i++) printf("%c",bdata_.head[i]); 

                     printf("\n");
                     printf("%s",iow_str);
                   }
                   eldata_.mct = 50;
                 } 

                 len = sprintf(iow_str,stress_format,
                   eldata_.n,eldata_.ma,sigav[0],sigav[1],sigav[2],sigav[3],sigav[4],sigav[5],
                   xx1,xx2,epsav[0],epsav[1],epsav[2],epsav[3],epsav[4],epsav[5]); 

                 fortiow_(iow_str,&len);

                 if(iofile_.ior < 0) printf("%s",iow_str);
                 
               }
               else 
                stcn13(ix,shps,xsj,sigl,&(_BLNK__.hr[prstrs_.nph-1]),
                                        &(_BLNK__.hr[prstrs_.nph+cdata_.numnp-1]),
                       lint,eldata_.nel,cdata_.numnp); 

           
              break;

          }

          default:
          break;
        }


}