/* $Header: /usr/home/sanjay/Feap/Elem/Ec/RCS/modl13.c,v 1.3 2000/01/28 20:33:06 sanjay Exp $ */

#include <math.h>

void modl13(double d[],
             double eps[6], double deps[6], double dt,
             double hold[6], double hnew[6],
             double sig[6], double ds[6][6]) {
/*
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      Simple viscoelastic model
c      with deviatoric relaxation in one exponential mode
c
c     Input:
c        d[0] - lam
c        d[1] - mu
c        d[2] - relaxation time
c        d[3] - viscosity strength
c        eps[6] - strains
c        deps[6] - strain increments
c        hold[6] - old history
c        dt      - time step
c
c     Output:
c        hnew[6]  - new history
c        sig[6]   - stresses
c        ds[6][6] - tangent
c-------------------------------------------------------------------
*/


      double lam,mu,tre,trs,p,muv,tau,k;
      double trde,xfact;
      int i,j;

      tre   = eps[0]+ eps[1]+ eps[2];
      trde = deps[0]+deps[1]+deps[2];

      lam    = d[0];
      mu     = d[1];
      tau    = d[2];
      muv    = d[3];

      xfact = exp(-dt/(2.0*tau));

      for(i=0;i<6;i++) {
       for(j=0;j<6;j++)   
        ds[i][j] = 0.0;
      }

      /* Elastic volumetric response */
      k = 2.0*mu/3.0 + lam;
      p = k*tre;
      
      sig[0] = p;
      sig[1] = p;
      sig[2] = p;
      sig[3] = 0.0;
      sig[4] = 0.0;
      sig[5] = 0.0;
 
      /* Viscous Deviatoric Response */

      for(i=0;i<3;i++) {
        eps[i]  =  eps[i] - tre/3.0;
        deps[i] = deps[i] - trde/3.0;
      }

      for(i=0;i<3;i++) {
         hnew[i] = xfact*xfact*hold[i] + xfact*deps[i];
         sig[i]  = sig[i]  + 2.0*mu*eps[i] + 2.0*muv*hnew[i];
      }
      for(i=3;i<6;i++) {
         hnew[i] = hold[i] + xfact*deps[i];
         sig[i]  = sig[i]  + mu*eps[i] + muv*hnew[i];
      }
     
      xfact = muv*xfact + mu;

      for(i=0; i<3; i++){
       for(j=0; j<3; j++) ds[i][j] = k-xfact*2.0/3.0;
      }

      ds[0][0] += 2.0*xfact;
      ds[1][1] += 2.0*xfact;
      ds[2][2] += 2.0*xfact;
      ds[3][3] += xfact;
      ds[4][4] += xfact;
      ds[5][5] += xfact;


}