elmt_shell_4n_q.c

利用语言编写的有限元分析软件

/* -----  Transform the loads and stiffness to global coords.  ----- */

       a = MatrixAllocIndirectDouble(3,3);
       b = dVectorAlloc(6);

       i0 = 0;
       for( ir=1; ir<=4; ir++ ) {
          for( ii=1; ii<=3; ii++ ) {
             b[ii-1] = 0.0;
             b[ii+2] = 0.0;
             for( k=1; k<=3; k++ ) {
                b[ii-1] += t[k-1][ii-1]*p[i0+k-1];
                b[ii+2] += t[k-1][ii-1]*p[i0+k+2];
             }
          }
          for( ii=1; ii<=6; ii++ ) {
             p[i0+ii-1] = b[ii-1];
          }
          j0 = i0;
          for( jc=ir; jc<=4; jc++ ) {
             i1 = i0;
             for( i=1; i<=2; i++ ) {
                j1 = j0;
                for( j=1; j<=2; j++ ) {
                   for( ii=1; ii<=3; ii++ ) {
                      for( jj=1; jj<=3; jj++ ) {
                         a[jj-1][ii-1] = 0.0;
                         for( k=1; k<=3; k++ ) {
                            a[jj-1][ii-1] += t[k-1][ii-1]*s[i1+k-1][jj+j1-1];
                         }
                      }
                   }
                   for( ii=1; ii<=3; ii++ ) {
                      for( jj=1; jj<=3; jj++ ) {
                         s[ii+i1-1][jj+j1-1] = 0.0;
                         for( k=1; k<=3; k++ ) {
                            s[ii+i1-1][jj+j1-1] += a[k-1][ii-1]*t[k-1][jj-1];
                         }
                      }
                   }
                j1 = j1 + 3;
                }
             i1 = i1 + 3;
             }

/* -----   Compute the symmetric block    -------------------------- */
 
             if( ir != jc ) {
                for( i=1; i<=6; i++ ) {
                   for( j=1; j<=6; j++ ) {
                      s[j0+j-1][i0+i-1] = s[i0+i-1][j0+j-1];
                   }
                }
             }
             j0 = j0 + ndf;
          }
          i0 = i0 + ndf;
       }

       MatrixFreeIndirectDouble(a,3);
       free ((char *) b);

#ifdef DEBUG
    printf("*** leaving rots06() ***\n");
#endif

}


/*===================================================================*/
/* Be transfered from rshp06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void rshp06( si, ss, tt, x, xsj, ndm )
double      **x;
double      ss, tt, *xsj;
int         si, ndm; 
{
/*  
double      s[4] = {-0.5, 0.5, 0.5, -0.5};
double      t[4] = {-0.5, -0.5, 0.5, 0.5};   
*/
double      **sx;
double      s2, t2, tp, xsj1;
int         i;
        
#ifdef DEBUG
    printf("*** Entering rshp06() ***\n");
#endif

/* -----  Form 4-node quatrilateral shape functions    ------------- */

       sx = MatrixAllocIndirectDouble( 2, 2 );

       shp[0][1][si] = 0.25*(1.0-tt);
       shp[0][2][si] = 0.25*(1.0+tt);
       shp[0][0][si] = -shp[0][1][si];
       shp[0][3][si] = -shp[0][2][si];
       shp[1][3][si] = 0.25*(1.0-ss);
       shp[1][2][si] = 0.25*(1.0+ss);
       shp[1][1][si] = -shp[1][2][si];
       shp[1][0][si] = -shp[1][3][si];
       shp[2][0][si] = shp[0][1][si]*(1.0-ss);
       shp[2][1][si] = shp[0][1][si]*(1.0+ss);
       shp[2][2][si] = shp[0][2][si]*(1.0+ss);
       shp[2][3][si] = shp[0][2][si]*(1.0-ss);
       s2        = (1.0-ss*ss)*0.5;
       t2        = (1.0-tt*tt)*0.5;
       shp[2][4][si] = s2*(1.0-tt);
       shp[2][5][si] = t2*(1.0+ss);
       shp[2][6][si] = s2*(1.0+tt);
       shp[2][7][si] = t2*(1.0-ss);
       shp[0][4][si] =-ss*(1.0-tt);
       shp[0][5][si] = t2;
       shp[0][6][si] =-ss*(1.0+tt);
       shp[0][7][si] =-t2;
       shp[1][4][si] =-s2;
       shp[1][5][si] =-tt*(1.0+ss);
       shp[1][6][si] = s2;
       shp[1][7][si] =-tt*(1.0-ss);

/* -----  Construct the jacobian and its inverse    ---------------- */

       sx[1][1] = x[0][0]*shp[0][0][si] + x[0][1]*shp[0][1][si] + 
                  x[0][2]*shp[0][2][si] + x[0][3]*shp[0][3][si];
       sx[0][1] = x[0][0]*shp[1][0][si] + x[0][1]*shp[1][1][si] + 
                  x[0][2]*shp[1][2][si] + x[0][3]*shp[1][3][si];
       sx[1][0] = x[1][0]*shp[0][0][si] + x[1][1]*shp[0][1][si] + 
                  x[1][2]*shp[0][2][si] + x[1][3]*shp[0][3][si];
       sx[0][0] = x[1][0]*shp[1][0][si] + x[1][1]*shp[1][1][si] + 
                  x[1][2]*shp[1][2][si] + x[1][3]*shp[1][3][si];

       *xsj     = sx[0][0]*sx[1][1] - sx[0][1]*sx[1][0];
       xsj1     = *xsj;

       if(xsj1 == 0.0 )
          xsj1 = 1.0;

       sx[1][1] = sx[1][1]/xsj1;
       sx[0][0] = sx[0][0]/xsj1;
       sx[0][1] =-sx[0][1]/xsj1;
       sx[1][0] =-sx[1][0]/xsj1;

/* -----  Form global derivatives    ------------------------------ */

       for( i=1; i<=8; i++ ) {
          tp              = shp[0][i-1][si]*sx[0][0] + shp[1][i-1][si]*sx[1][0];
          shp[1][i-1][si] = shp[0][i-1][si]*sx[0][1] + shp[1][i-1][si]*sx[1][1];
          shp[0][i-1][si] = tp; 
       }

/* -----  Form the rotational and 5th shape functions  ------------- */

       hshp06 ( si );

        MatrixFreeIndirectDouble(sx,2);

#ifdef DEBUG
    printf("*** leaving rshp06() ***\n");
#endif

}


/*===================================================================*/
/* Be transfered from stre06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void stre06( p, d, vl, ndm, nel, k, xx, yy, zz, eps, sigi, sigb )
ARRAY       *p;
double      *d, **vl, *eps, *sigi, *sigb;
double      *xx, *yy, *zz;
int         ndm, nel, k; 
{
int         i, j;

#ifdef DEBUG
    printf("*** Entering stre06() ***\n");
#endif

/*
#ifdef DEBUG
    printf("d vector\n");
    for(i=1; i<=20; i++) printf("%5d%15.8e\n", i, d[i-1]);
    printf("p->coord array\n");
           printf("                   1              2              3               4\n");
    for(i=1; i<=3; i++){ 
        printf("%5d", i);
        for(j=1; j<=4; j++) 
           printf("%15.4e", p->coord[i-1][j-1].value);
        printf("\n");
    }

    dMatrixPrint("vl", vl, 6, 4);
#endif
*/

/* -----  Compute the membrane and bending strains    -------------- */

      dktq06( k-1 );

      for( i=1; i<=6; i++ ) 
         eps[i-1] = 0.0;

      *xx = 0.0;
      *yy = 0.0;
      *zz = 0.0;

      for( j=1; j<=nel; j++ ) {

         *xx += shp[2][j-1][k-1]*p->coord[0][j-1].value;
         *yy += shp[2][j-1][k-1]*p->coord[1][j-1].value;
         *zz += shp[2][j-1][k-1]*p->coord[2][j-1].value;
         eps[0] += shp[0][j-1][k-1]*vl[0][j-1] - shp1[0][j-1][k-1]*vl[5][j-1];
         eps[2] += shp[1][j-1][k-1]*vl[1][j-1] - shp2[1][j-1][k-1]*vl[5][j-1];
         eps[1] += shp[0][j-1][k-1]*vl[1][j-1] + shp[1][j-1][k-1]*vl[0][j-1] 
                   - (shp1[1][j-1][k-1] + shp2[0][j-1][k-1])*vl[5][j-1];

         for( i=ii1; i<=ii2; i++ ) {
            eps[3] += bm[0][i-1][j-1]*vl[i-1][j-1];
            eps[5] += bm[1][i-1][j-1]*vl[i-1][j-1];
            eps[4] += bm[2][i-1][j-1]*vl[i-1][j-1];
         }
      }

      sigi[0] = d[0]*eps[0] + d[1]*eps[2];
      sigi[2] = d[0]*eps[2] + d[1]*eps[0];
      sigi[1] = d[2]*eps[1];

      sigi = pstres06(sigi);

      sigb[0] = d[3]*eps[3] + d[4]*eps[5];
      sigb[2] = d[4]*eps[3] + d[3]*eps[5];
      sigb[1] = d[5]*eps[4];

      sigb = pstres06(sigb);

#ifdef DEBUG
    printf("*** leaving stre06() ***\n");
#endif

}
      

/* ----------------------------------------*/
/*     shp_prt                             */
/* ----------------------------------------*/

void shp_prt(nel)
int nel;
{
int i, j, k;

    printf("shp[3][8][9] array\n");
    for(k=1; k<=3; k++) {
    printf("k = %2d\n           1         2         3         4         5         6         7         8\n", k);
       for(i=1; i<=9; i++){ 
          printf("%2d", i);
          for(j=1; j<=8; j++) 
             printf("%10.3e", shp[k-1][j-1][i-1]);
          printf("\n");
       }
    }
    printf("shp1[3][4][9] array\n");
    for(k=1; k<=3; k++) {
    printf("k = %2d\n            1          2          3          4\n", k);
       for(i=1; i<=9; i++){ 
          printf("%2d", i);
          for(j=1; j<=4; j++) 
             printf("%11.3e", shp1[k-1][j-1][i-1]);
          printf("\n");
       }
    }
    printf("shp2[3][4][9] array\n");
    for(k=1; k<=3; k++) {
    printf("k = %2d\n            1          2          3          4\n", k);
       for(i=1; i<=9; i++){ 
          printf("%2d", i);
          for(j=1; j<=4; j++) 
             printf("%11.3e", shp2[k-1][j-1][i-1]);
          printf("\n");
       }
    }
    printf("bm[3][6][9] array\n");
    for(k=1; k<=3; k++) {
    printf("k = %2d\n            1          2          3          4          5          6\n", k);
       for(i=1; i<=9; i++){ 
          printf("%2d", i);
          for(j=1; j<=6; j++) 
             printf("%11.3e", bm[k-1][j-1][i-1]);
          printf("\n");
       }
    }

    printf("nel = %3d, ii1 = %3d, ii2 = %3d\n",nel,ii1,ii2);

}


/*===================================================================*/
/* Calculation  of  Transformation  Array  and  Surface  Coordiantes */
/* Be transfered from tran06 by L. Jin,  Jan 14                      */
/* ==================================================================*/
     
void tran06(p, yl, t)
ARRAY      *p;
double     **yl, **t;
{
int      i, j, k;
double   v1, v2, htol, d11, d12, *x0;

#ifdef DEBUG
    printf("*** Entering tran06() ***\n");
#endif

/* --- Compute the inplane direction cosines (bisect diagonals) --- */

       x0 = dVectorAlloc(3);

       for( i=1; i<=3; i++ ) {
           t[0][i-1] = p->coord[i-1][2].value - p->coord[i-1][0].value;
           t[1][i-1] = p->coord[i-1][1].value - p->coord[i-1][3].value;
       }

       d11 = sqrt(t[0][0]*t[0][0] + t[0][1]*t[0][1] + t[0][2]*t[0][2]);
       d12 = sqrt(t[1][0]*t[1][0] + t[1][1]*t[1][1] + t[1][2]*t[1][2]);
                   
       for( i=1; i<=3; i++ ) {
           v1 = t[0][i-1]/d11;
           v2 = t[1][i-1]/d12;
           t[0][i-1] = v1 + v2;
           t[1][i-1] = v1 - v2;
       }

       d11 = sqrt(t[0][0]*t[0][0] + t[0][1]*t[0][1] + t[0][2]*t[0][2]);
       d12 = sqrt(t[1][0]*t[1][0] + t[1][1]*t[1][1] + t[1][2]*t[1][2]);

       for( i=1; i<=3; i++ ) {
           t[0][i-1] = t[0][i-1]/d11;
           t[1][i-1] = t[1][i-1]/d12;

/* ------------ Compute the center (0,0) displacement ------------- */

           x0[i-1] = 0.25*(p->coord[i-1][0].value + p->coord[i-1][1].value +
                           p->coord[i-1][2].value + p->coord[i-1][3].value);
       }

/* -------------- Compute the normal to the surface --------------- */

       t[2][0] = t[0][1]*t[1][2] - t[1][1]*t[0][2] ;
       t[2][1] = t[0][2]*t[1][0] - t[1][2]*t[0][0] ;
       t[2][2] = t[0][0]*t[1][1] - t[1][0]*t[0][1] ;

/* ------- Compute the projected middle surface coordinates ------- */

       for( i=1; i<=4; i++ ) {
          for( j=1; j<=3; j++ ) {   
             yl[j-1][i-1] = 0.0;
             for( k=1; k<=3; k++ ) {
                yl[j-1][i-1] += t[j-1][k-1]*(p->coord[k-1][i-1].value - x0[k-1]);
             }
          }
       }

/* - Set offset coordinates to zero if small compared to plan size - */

       htol = 0.0;

       for( i=1; i<=4; i++ ) {
          htol = MAX( htol, ABS(yl[0][i-1]));
          htol = MAX( htol, ABS(yl[1][i-1]));
       }

       htol = htol * 1.0e-7;

       for( i=1; i<=4; i++ ) {
          if( ABS(yl[2][i-1]) <= htol )
             yl[2][i-1] = 0.0;
       }

       free ((char *) x0);

#ifdef DEBUG
    printf("*** leaving tran06() ***\n");
#endif

}


/* ============================================= */
/*          pstres06                             */
/* ============================================= */

double *pstres06(sig)
double *sig;
{
double xi1, xi2, rho;

#ifdef DEBUG
     printf("*** entering pstres06() ***\n");
#endif

     xi1 = 0.5*(sig[0] + sig[2]);
     xi2 = 0.5*(sig[0] - sig[2]);
     rho = sqrt(xi2*xi2 + sig[1]*sig[1]);
     sig[3] = xi1 + rho;
     sig[4] = xi1 - rho;
     sig[5] = 45.0;
     if(xi2 != 0.0) sig[5] = 22.5*atan2(sig[1],xi2)/atan(1.0);
     
#ifdef DEBUG
    printf("*** leaving pstres06() ***\n");
#endif

     return (sig);

}