#include <stdio.h>
/* define s(nst,nst) nad xl(ndm,nel) as a zero
index arrays (column major) based
upon pointers that come
from the FORTRAN parts of FEAP*/
#define s(a,b) (*(sf+(a)+((b)*nst)))
#define xl(a,b) (*(xlf+(a)+((b)*ndm)))
struct {
double hr[1];
int mr[1000];
} _BLNK__; /* The blank common structure
*/
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 string for stress output header */
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 string for stress output formats */
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 your typical fortran elements
*/
double shps[4][4][3]; /*
swapped order from your typical fortran elements*/
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:
/* Integration points and weights */
tmp1 =
2; /* Gauss points per direction */
lint =
4; /* Gauss points in element
*/
int2d_(&tmp1,&lint,*sg);
/* Uses tmp1 and tmp2 for pointer extraction
since this is a call to a FORTRAN routine
Note the last argument is *sg which is actually a pointer to a double precision
number. */
/*
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. *(sg+l) is a pointer
to the l'th Gauss point's coordinates and weight.
c
**(shps+l) is a pointer to the l'th Gauss points shape function derivatives
and value. *(eps+l) is a pointer to the
c
strains for the l'th Gauss point.
*/
for(l=0;
l<lint; l++)
shp2d_(*(sg+l), xlf, **(shps+l), xsj+l,&ndm,&eldata_.nel,ix,&FALSE);
/* To pass a false to FORTRAN use a pointer to a
zero */
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++) {
/* Note the shift on the arguments into the
blank common arrays to account for zero based indexing in C */
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;
}
/* Compute internal forces */
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;
}
/* Do tangent for isw = 3 */
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;
/* Stress outputs */
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) {
/* Use fortiow_() cover function for outputs */
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);
}
/* Stress plots */
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;
}
}