elmt_shell_4n_q.c

有限元分析源代码

#ifdef DEBUG
    printf(" In elmt_shell_4nodes_q(): \n");
    printf("    : leaving case of MASS_MATRIX\n");
#endif
        
        break;
           
        default:
        break;
    }

/*
        MatrixFreeIndirectDouble(yl, 3);
        MatrixFreeIndirectDouble(tr, 3);
        MatrixFreeIndirectDouble(gshp1, 3);
        MatrixFreeIndirectDouble(gshp2, 3);
        MatrixFreeIndirectDouble(vl, 6);
        MatrixFreeIndirectDouble(btd, 3);
        MatrixFreeIndirectDouble(s, nst);
 
        free ((char *) eps);
        free ((char *) sigi);
        free ((char *) sigb);
        free ((char *) sg);
        free ((char *) tg);
        free ((char *) wg);
        free ((char *) gg);
        free ((char *) dvl);
        free ((char *) pp);
        free ((char *) d);
        free ((char *) bl);
*/

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

    return(p);
}


/*===================================================================*/
/* Be transfered from dktb06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void dktb06( sg, tg,xsj )
double sg, tg,xsj;
{
double      *el, **shn, *bd, *cd;
double      shm1, shm2;
double      *a1, *b1, *c1, *d1, *e1;
double      *a2, *b2, *c2, *d2, *e2;
int         i, j, k;
  
#ifdef DEBUG
    printf("*** Entering dktb06() ***\n");
#endif

       el  = dVectorAlloc(3);
       bd  = dVectorAlloc(3);
       cd  = dVectorAlloc(3);
       a1  = dVectorAlloc(3);
       a2  = dVectorAlloc(3);
       b1  = dVectorAlloc(3);
       b2  = dVectorAlloc(3);
       c1  = dVectorAlloc(3);
       c2  = dVectorAlloc(3);
       d1  = dVectorAlloc(3);
       d2  = dVectorAlloc(3);
       e1  = dVectorAlloc(3);
       e2  = dVectorAlloc(3);
       shn = MatrixAllocIndirectDouble( 2, 3 );

/* -------    Compute the area coordinates for the point    -------- */

       el[0] = 0.25 * (1.0 - sg) * (1.0 - tg);
       el[1] = 0.25 * (1.0 + sg) * (1.0 - tg);
       el[2] = 0.5  * (1.0 + tg);

/* -------     Form the shape function terms for trains     -------- */

       for( i=1; i<=3; i++ ) {
          bd[i-1] = b[i-1]/xsj;
          cd[i-1] = c[i-1]/xsj;
       }
       for( i=1; i<=3; i++ ) {
          j = i%3 + 1;
          k = j%3 + 1;
          shn[0][i-1] = bd[i-1] * (el[i-1] - 1.0);
          shn[1][i-1] = cd[i-1] * (el[i-1] - 1.0);
          shm1 = bd[j-1]*el[k-1] + bd[k-1]*el[j-1];
          shm2 = cd[j-1]*el[k-1] + cd[k-1]*el[j-1];
          a1[i-1] = aa[i-1] * shm1;
          a2[i-1] = aa[i-1] * shm2;
          b1[i-1] = bb[i-1] * shm1;
          b2[i-1] = bb[i-1] * shm2;
          c1[i-1] = cc[i-1] * shm1;
          c2[i-1] = cc[i-1] * shm2;
          d1[i-1] = dd[i-1] * shm1;
          d2[i-1] = dd[i-1] * shm2;
          e1[i-1] = ee[i-1] * shm1;
          e2[i-1] = ee[i-1] * shm2;
       }

/* ---------      Plate train displacement matrix       ------------ */

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

          j = i%3 + 1;
          k = j%3 + 1;
          bm[0][2][i-1] = a1[k-1] - a1[j-1];
          bm[0][3][i-1] = b1[k-1] + b1[j-1];
          bm[0][4][i-1] = c1[k-1] + c1[j-1] - shn[0][i-1];
          bm[1][2][i-1] = d2[k-1] - d2[j-1];
          bm[1][3][i-1] =-e2[k-1] + e2[j-1] + shn[1][i-1];
          bm[1][4][i-1] =-b2[k-1] + b2[j-1];
          bm[2][2][i-1] = a2[k-1] - a2[j-1] + d1[k-1] - d1[j-1];
          bm[2][3][i-1] =-e1[k-1] - e1[j-1] + shn[0][i-1] - bm[1][4][i-1];
          bm[2][4][i-1] = c2[k-1] + c2[j-1] - shn[1][i-1] - bm[0][3][i-1];
       }
    
       free ((char *) el);
       free ((char *) bd);
       free ((char *) cd);
       MatrixFreeIndirectDouble(shn,2);
       free ((char *) a1);
       free ((char *) a2);
       free ((char *) b1);
       free ((char *) b2);
       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

}


/*===================================================================*/
/* Be transfered from dktq06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void dktq06( si )
int         si;
{
int         i, j, k;
  
#ifdef DEBUG
    printf("*** Entering dktq06() ***\n");
#endif

    /* 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

}


/*===================================================================*/
/* Be transfered from hshp06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void hshp06( si )
int         si;
{
double      *shx, *shy;
int         i, j, k;
  
#ifdef DEBUG
    printf("*** Entering hshp06() ***\n");
#endif

/* -----    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;
          }
       }
    
       free ((char *) shx );
       free ((char *) shy );
       
#ifdef DEBUG
    printf("*** leaving hshp06() ***\n");
#endif

}


/*===================================================================*/
/* Be transfered from jacq06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
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

/* -----    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

}


/*===================================================================*/
/* Be transfered from jtri06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void jtri06( xl, xsj )
double      **xl, *xsj;
{
double      sql, cs, bs;
int         i, j, k;
  
#ifdef DEBUG
    printf("*** Entering jtri06() ***\n");
#endif

/* -----    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;

/* -----   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

}


/*===================================================================*/
/* Be transfered from proj06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
void proj06( i1, j1, s, zi, zj, nst )
double      **s;
double      zi, zj;
int         i1, j1, nst; 
{
int         i;
  
#ifdef DEBUG
    printf("*** Entering proj06() ***\n");
#endif

/* -----    Modify the stiffness for offset projections      ------- */
/* -----    Postmultiply 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  ];
       }

/* --- Premultiply using modified terms from prosmultiplication  --- */

       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];
       }

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

}


/*===================================================================*/
/* Be transfered from rots06 by L. Jin,  Jan 16                      */
/* ==================================================================*/
     
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