elmt_shell_4n_q.c

有限元程序

   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);

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

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

   /* [d] : Form the rotational and 5th shape functions */

   hshp06 ( si );

   MatrixFreeIndirectDouble(sx,2);

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

}

/*
 *  ========================================================
 *  stre06() : Compute the membrane and bending strains
 *  
 *  Input : 
 *  Output : 
 *  ========================================================
 */
     
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

    /* [a] : 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() : print shape function arrays.
 *  
 *  Input  : 
 *  Output : void.
 *  ======================================================================
 */

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

    /* [a] : print shp[3][8][9] array */

    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");
       }
    }

    /* [b] : print shp1[3][4][9] array */

    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");
       }
    }

    /* [c] : print shp2[3][4][9] array */

    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");
       }
    }

    /* [d] : print bm[3][6][9] array */

    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);

}


/*
 *  ======================================================================
 *  tran06() : Calculation of Transformation Array and Surface Coordinates
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */
     
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

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

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

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

    /* [d] : 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() : 
 * 
 *  Input  : 
 *  Output : 
 *  ======================================================
 */ 

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

/*
 *  =====================================================================
 *  print_property_shell_4n_q() : Print SHELL_4NQ element properties
 *
 *  Input  : 
 *  Output : 
 *  =====================================================================
 */

#ifdef __STDC__
void print_property_shell_4n_q(EFRAME *frp, int i)
#else
void print_property_shell_4n_q(frp, i)
EFRAME    *frp;
int          i;                 /* elmt_attr_no */
#endif
{
int     UNITS_SWITCH;
ELEMENT_ATTR    *eap;

     UNITS_SWITCH = CheckUnits();
     eap = &frp->eattr[i-1];

     if( PRINT_MAP_DOF == ON ) {
        if(frp->no_dof == 3 || frp->no_dof == 2) { 
           printf("             ");
           printf("         : gdof [0] = %4d : gdof[1] = %4d : gdof[2] = %4d\n",
                           eap->map_ldof_to_gdof[0],
                           eap->map_ldof_to_gdof[1],
                           eap->map_ldof_to_gdof[2]);
        }

        if(frp->no_dof == 6) { /* 3d analysis */
           printf("             ");
           printf("         : dof-mapping : gdof[0] = %4d : gdof[1] = %4d : gdof[2] = %4d\n",
                           eap->map_ldof_to_gdof[0],
                           eap->map_ldof_to_gdof[1],
                           eap->map_ldof_to_gdof[2]);
           printf("             ");
           printf("                         gdof[3] = %4d : gdof[4] = %4d : gdof[5] = %4d\n",
                           eap->map_ldof_to_gdof[3],
                           eap->map_ldof_to_gdof[4],
                           eap->map_ldof_to_gdof[5]);
        } 
     }

     switch(UNITS_SWITCH) {
       case ON:
        UnitsSimplify( eap->work_material[0].dimen );
        UnitsSimplify( eap->work_material[2].dimen );
        UnitsSimplify( eap->work_material[5].dimen );
        UnitsSimplify( eap->work_section[2].dimen );
        UnitsSimplify( eap->work_section[10].dimen );
        if( eap->work_material[0].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Young's Modulus =  E = %16.3e %s\n",
                           eap->work_material[0].value/eap->work_material[0].dimen->scale_factor,
                           eap->work_material[0].dimen->units_name);
        }
        if( eap->work_material[4].value != 0.0 ) {
           printf("             ");
           printf("         : Poisson's ratio = nu = %16.3e   \n", eap->work_material[4].value);
        }
        if( eap->work_material[2].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Yielding Stress = fy = %16.3e %s\n",
                           eap->work_material[2].value/eap->work_material[2].dimen->scale_factor,
                           eap->work_material[2].dimen->units_name);
        }
	if( eap->work_material[5].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Density         = %16.3e %s\n",
                           eap->work_material[5].value/eap->work_material[5].dimen->scale_factor,
                           eap->work_material[5].dimen->units_name);
	}
	if( eap->work_section[2].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Inertia Izz     = %16.3e %s\n",
                           eap->work_section[2].value/eap->work_section[2].dimen->scale_factor,
                           eap->work_section[2].dimen->units_name);
	}
	if( eap->work_section[10].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Area            = %16.3e %s\n",
                           eap->work_section[10].value/eap->work_section[10].dimen->scale_factor,
                           eap->work_section[10].dimen->units_name);
	}
       break;
       case OFF:
         printf("             ");
         printf("         : Young's Modulus =  E = %16.3e\n",
                            eap->work_material[0].value);
         printf("             ");
         printf("         : Yielding Stress = fy = %16.3e\n",
                            eap->work_material[2].value);
         printf("             ");
         printf("         : Poisson's ratio = nu = %16.3e   \n", eap->work_material[4].value);
         printf("             ");
         printf("         : Density         = %16.3e\n",
                            eap->work_material[5].value);
         printf("             ");
         printf("         : Inertia Izz     = %16.3e\n",
                            eap->work_section[2].value);
         printf("             ");
         printf("         : Area            = %16.3e\n",
                            eap->work_section[10].value);
        break;
        default:
        break;
     }
}