elmt_shell_4n_q.c

有限元程序

   free ((char *) c1);
   free ((char *) c2);
   free ((char *) d1);
   free ((char *) d2);
   free ((char *) e1);
   free ((char *) e2);

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

}


/*
 *  ======================================================================
 *  dktq06() : Form strain-displacement array for DKQ element
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */
     
void dktq06( si )
int         si;
{
int         i, j, k;
  
#ifdef DEBUG
       printf("*** Entering dktq06() ***\n");
#endif

   /* [a] : Form strain-displacement array for DKQ element */

   for( i=1; i<=4; i++ ) {

       j = (i+2) % 4 + 1;
       bm[0][2][i-1] = aa[i-1]*shp[0][i+3][si] - aa[j-1]*shp[0][j+3][si];
       bm[0][3][i-1] = bb[i-1]*shp[0][i+3][si] + bb[j-1]*shp[0][j+3][si];
       bm[0][4][i-1] = cc[i-1]*shp[0][i+3][si] + cc[j-1]*shp[0][j+3][si] - shp[0][i-1][si];
       bm[1][2][i-1] = dd[i-1]*shp[1][i+3][si] - dd[j-1]*shp[1][j+3][si];
       bm[1][3][i-1] =-ee[i-1]*shp[1][i+3][si] - ee[j-1]*shp[1][j+3][si] + shp[1][i-1][si];
       bm[1][4][i-1] =-bb[i-1]*shp[1][i+3][si] - bb[j-1]*shp[1][j+3][si];
       bm[2][2][i-1] = aa[i-1]*shp[1][i+3][si] - aa[j-1]*shp[1][j+3][si] +
                       dd[i-1]*shp[0][i+3][si] - dd[j-1]*shp[0][j+3][si];
       bm[2][3][i-1] =-ee[i-1]*shp[0][i+3][si] - ee[j-1]*shp[0][j+3][si] +
                               shp[0][i-1][si] - bm[1][4][i-1];
       bm[2][4][i-1] = cc[i-1]*shp[1][i+3][si] + cc[j-1]*shp[1][j+3][si] -
                               shp[1][i-1][si] - bm[0][3][i-1];
   }
    
#ifdef DEBUG
       printf("*** leaving dktq06() ***\n");
#endif

}


/*
 *  ======================================================================
 *  hshp06() : Form the shape function 
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

void hshp06( si )
int si;
{
double *shx, *shy;
int    i, j, k;
  
#ifdef DEBUG
       printf("*** Entering hshp06() ***\n");
#endif

   /* [a] : Form the shape function */

   shx = dVectorAlloc( 4 );
   shy = dVectorAlloc( 4 );

   for( k=1; k<=3; k++ ) {
        for( i=1; i<=4; i++ ) {
             shx[i-1] = shp[k-1][i+3][si]*c[i-1]; 
             shy[i-1] = shp[k-1][i+3][si]*b[i-1]; 
        }

        j = 4;
        for( i=1; i<=4; i++ ) {
             shp1[k-1][i-1][si] = shy[i-1] - shy[j-1];
             shp2[k-1][i-1][si] = shx[i-1] - shx[j-1];
             j = i;
        }
   }

   /* [b] : free working memory */
    
   free ((char *) shx );
   free ((char *) shy );
       
#ifdef DEBUG
       printf("*** leaving hshp06() ***\n");
#endif

}


/*
 *  ======================================================================
 *  jacq06() : Form geometric constants for DKQ element.
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

void jacq06( xl )
double **xl;
{
double  cxx, cyy, cxy, dz, ad, a4, b2, c2, sql, x31, y31, x42, y42;
int     i, j, k;
  
#ifdef DEBUG
       printf("*** Entering jacq06() ***\n");
#endif

    /* [a] : Form geometric constants for DKQ element   ------- */

    cxx = 0.0;
    cyy = 0.0;     
    cxy = 0.0;     

    for( i=1; i<=4; i++ ) {
       k = i % 4 + 1;
       b[i-1] = xl[1][k-1] - xl[1][i-1];
       c[i-1] = xl[0][i-1] - xl[0][k-1];
       dz     = xl[2][k-1] - xl[2][i-1];
       b2     = b[i-1]*b[i-1];
       c2     = c[i-1]*c[i-1];
       cxx   += dz*b2/(b2 + c2);
       cyy   += dz*c2/(b2 + c2);
       cxy   += (dz+dz)*b[i-1]*c[i-1]/(b2 + c2);
       sql    = (b2 + c2)/0.75;
       aa[i-1] = (c[i-1] + c[i-1])/sql;
       bb[i-1] =  b[i-1] * c[i-1] /sql;
       cc[i-1] =  c2/sql;
       dd[i-1] =-(b[i-1] + b[i-1])/sql;
       ee[i-1] =  b2/sql;
        b[i-1] =  b[i-1]/8.0;
        c[i-1] =  c[i-1]/8.0;
    }

    if( xl[2][0] != 0.0 ) {
        x31 = xl[0][2] - xl[0][0];
        y31 = xl[1][2] - xl[1][0];
        x42 = xl[0][3] - xl[0][1];
        y42 = xl[1][3] - xl[1][1];
        a4  = (x31*y42 - x42*y31)*2.0;
        ad  = a4*a4/4.0;
        x31 = x31/ad;
        y31 = y31/ad;
        x42 = x42/ad;
        y42 = y42/ad;

        bm[0][0][0] = -cxx*x42;
        bm[1][0][0] = -cyy*x42;
        bm[2][0][0] = -cxy*x42;

        bm[0][1][0] = -cxx*y42;
        bm[1][1][0] = -cyy*y42;
        bm[2][1][0] = -cxy*y42;

        bm[0][5][0] = -cxx/a4;
        bm[1][5][0] = -cyy/a4;
        bm[2][5][0] = -cxy/a4;

        bm[0][0][1] =  cxx*x31;
        bm[1][0][1] =  cyy*x31;
        bm[2][0][1] =  cxy*x31;

        bm[0][1][1] =  cxx*y31;
        bm[1][1][1] =  cyy*y31;
        bm[2][1][1] =  cxy*y31;

        bm[0][5][1] = -cxx/a4;
        bm[1][5][1] = -cyy/a4;
        bm[2][5][1] = -cxy/a4;

        for( i=1; i<=3; i++ ) {
           bm[i-1][0][2] = -bm[i-1][0][0];
           bm[i-1][1][2] = -bm[i-1][1][0];
           bm[i-1][5][2] =  bm[i-1][5][0];

           bm[i-1][0][3] = -bm[i-1][0][1];
           bm[i-1][1][3] = -bm[i-1][1][1];
           bm[i-1][5][3] =  bm[i-1][5][1];
        }
    }
    
#ifdef DEBUG
       printf("*** Leaving jacq06() ***\n");
#endif

}

/*
 *  ======================================================================
 *  jtri06() : Form the terms for Jacobian of a triangle 
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

void jtri06( xl, xsj )
double      **xl, *xsj;
{
double      sql, cs, bs;
int         i, j, k;
  
#ifdef DEBUG
       printf("*** Entering jtri06() ***\n");
#endif

    /* [a] : Form the terms for the Jacobian of a triangle */

    for( i=1; i<=3; i++ ) {
        j = i%3 + 1;
        k = j%3 + 1;
        b[i-1] = xl[1][k-1] - xl[1][j-1];
        c[i-1] = xl[0][j-1] - xl[0][k-1];
        sql    = 4.0*(b[i-1]*b[i-1] + c[i-1]*c[i-1]);
        cs     = c[i-1]/sql;
        bs     = b[i-1]/sql;
        aa[i-1] = 6.0*cs;
        bb[i-1] = 3.0*bs*c[i-1];
        cc[i-1] = c[i-1]*cs - b[i-1]*(bs+bs);
        dd[i-1] =-6.0*bs;
        ee[i-1] = b[i-1]*bs - c[i-1]*(cs+cs);
    }

    b[3]  = 0.0;
    c[3]  = 0.0;
    aa[3] = 0.0;
    bb[3] = 0.0;
    cc[3] = 0.0;
    dd[3] = 0.0;
    ee[3] = 0.0;

    /* [b] : Store jacobian */

    *xsj   = - xl[0][0]*b[0] - xl[0][1]*b[1] - xl[0][2]*b[2];

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

}

/*
 *  ======================================================================
 *  proj06() : Modify stiffness for offset projections 
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

void proj06( i1, j1, s, zi, zj, nst )
double      **s;
double      zi, zj;
int         i1, j1, nst; 
{
int         i;
  
   /* [a] : Post-multiply by transformation */

   for( i=1; i<=6; i++ ) {
        s[i1+i-1][j1+3] += zj*s[i1+i-1][j1+1];
        s[i1+i-1][j1+4] -= zj*s[i1+i-1][j1  ];
   }

   /* [b] : Premultiply using modified terms from post multiplication */

   for( i=1; i<=6; i++ ) {
        s[i1+3][j1+i-1] += zi*s[i1+1][j1+i-1];
        s[i1+4][j1+i-1] -= zi*s[i1  ][j1+i-1];
   }
}

/*
 *  ======================================================================
 *  rots06() : Transform the loads and stiffness to global coords.
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

void rots06( s, p, t, nst, ndf )
double      **s, *p, **t;
int         nst, ndf; 
{
double **a, *b;
int    i, i0, i1, ir, ii, j, j0, j1, jc, jj, k;

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

   /* [a] : 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;
   }

   /* [b] : Free memory for working matrices */

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

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

}

/*
 *  ======================================================================
 *  rshp06() : Form 4-node quatrilateral shape functions
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */
     
void rshp06( si, ss, tt, x, xsj, ndm )
double  **x;
double  ss, tt, *xsj;
int     si, ndm; 
{
double  **sx;
double  s2, t2, tp, xsj1;
int     i;
        
#ifdef DEBUG
       printf("*** Entering rshp06() ***\n");
#endif

   /* [a] : 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);