elmt_plate.c

有限元程序

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

            H_Print = H_Print-1;
            if(H_Print <= 0) {
               H_Print = 50;
               printf("---------------------------------------------------------------------------\n");
               printf("DKQ Plate Bending Element\n");
               printf("elmt   mat   x-coord   z-coord      Mxx/length     Mxz/length    Mzz/legnth\n");
               printf("---------------------------------------------------------------------------\n");
               if(UNITS_SWITCH == ON) {   
                  printf("Units");
	          switch( UnitsType ) {
	            case US:
                       dimen = DefaultUnits("lbf");
                       break;
	            case SI:
	            case SI_US:
	            default:
                       dimen = DefaultUnits("N");
                       break;
                  }
                  printf("           %s        %s           %s            %s             %s",
                         p->coord[0][0].dimen->units_name,
                         p->coord[2][0].dimen->units_name,
                         dimen->units_name, dimen->units_name, dimen->units_name);
                  printf("\n");
               }
            }

            printf("%4d %s %9.3f %9.3f", p->elmt_no, p->material_name,
                   xx/p->coord[0][0].dimen->scale_factor,
                   yy/p->coord[2][0].dimen->scale_factor);

            printf("%14.5e %14.5e %14.5e\n",
                   sig[0]/dimen->scale_factor,
                   sig[1]/dimen->scale_factor,
                   sig[2]/dimen->scale_factor);
            printf("\n");
            free((char *) dimen->units_name);
            free((char *) dimen);

            break;
       case LOAD_MATRIX:   /* Compute The Element Residual Vector */
            break;
    }

    if(isw != PROPTY){
       MatrixFreeIndirectDouble(shp, 3);
       MatrixFreeIndirectDouble(bm, 3);
    }

    return(p);
}


/* ============= */
/* Jacqud(x,ndm) */
/* ============= */

int jacqud(x,aa,bb,cc,dd,ee)
QUANTITY **x;
double aa[5],bb[5],cc[5],dd[5],ee[5];
{
int i,k;
double b,c,sql;

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

       b = x[2][k-1].value-x[2][i-1].value;
       c = x[0][i-1].value-x[0][k-1].value;
       sql = b*b+c*c;

       aa[i-1] = 1.5*c/sql;
       bb[i-1] = 0.75*b*c/sql;
       cc[i-1] = (0.25*c*c-0.5*b*b)/sql;
       dd[i-1] = -1.5*b/sql;
       ee[i-1] = (0.25*b*b-0.5*c*c)/sql;
   }

   return;
}


/* ============================= */
/* dktqbm(shm,bm,aa,bb,cc,dd,ee) */
/* ============================= */

double **dktqbm(shm,bm,aa,bb,cc,dd,ee)
double **shm,**bm;
double aa[5],bb[5],cc[5],dd[5],ee[5];
{
int i,j,i1,i2,i3;

   i1 = 1;
   for(i=1; i<=4; i++) {
      j  = (i+2)%4 +1;
      i2 = i1+1;
      i3 = i2+1;
      
      bm[0][i1-1] =  aa[i-1]*shm[0][i+4-1] - aa[j-1]*shm[0][j+4-1];
      bm[0][i2-1] =  bb[i-1]*shm[0][i+4-1] + bb[j-1]*shm[0][j+4-1];
      bm[0][i3-1] =  cc[i-1]*shm[0][i+4-1] + cc[j-1]*shm[0][j+4-1] - shm[0][i-1];
      bm[1][i1-1] =  dd[i-1]*shm[1][i+4-1] - dd[j-1]*shm[1][j+4-1];
      bm[1][i2-1] = -ee[i-1]*shm[1][i+4-1] - ee[j-1]*shm[1][j+4-1] + shm[1][i-1];
      bm[1][i3-1] = -bb[i-1]*shm[1][i+4-1] - bb[j-1]*shm[1][j+4-1];
      bm[2][i1-1] =  aa[i-1]*shm[1][i+4-1] - aa[j-1]*shm[1][j+4-1]
                    +dd[i-1]*shm[0][i+4-1] - dd[j-1]*shm[0][j+4-1];
      bm[2][i2-1] = -ee[i-1]*shm[0][i+4-1] - ee[j-1]*shm[0][j+4-1] + shm[0][i-1] - bm[1][i3-1];
      bm[2][i3-1] =  cc[i-1]*shm[1][i+4-1] + cc[j-1]*shm[1][j+4-1] - shm[1][i-1] - bm[0][i2-1];
      i1 = i1+3;
   }

   return(bm);
} 


/* ========================== */
/* qushp8(s,t,shp,x,xsj)      */
/* ========================== */

double **qushp8(s,t,shp,x,xsj) 
double           s,t;
double         **shp;
QUANTITY    **x;
double          *xsj;
{
int i;
double xs,xt,ys,yt,ss,tt,sn,tn,si[5],ti[5]; 

#ifdef DEBUG
       printf("*** In qushp8() : s   = %10.5f t   = %10.5f\n", s,t);
#endif

   si[0] = -1.0;
   si[1] =  1.0;
   si[2] =  1.0;
   si[3] = -1.0;

   ti[0] = -1.0;
   ti[1] = -1.0;
   ti[2] =  1.0;
   ti[3] =  1.0;

   xs = 0.0; xt = 0.0;
   ys = 0.0; yt = 0.0;

   for(i=1;i<=4;i++) {
       ss = si[i-1]*s;
       tt = ti[i-1]*t;
       sn = si[i-1]*(1.+tt);
       tn = ti[i-1]*(1.+ss);
       xs = xs+sn*x[0][i-1].value;
       xt = xt+tn*x[0][i-1].value;
       ys = ys+sn*x[2][i-1].value;
       yt = yt+tn*x[2][i-1].value;

       shp[0][i-1] = 0.25*sn*(ss+ss+tt);
       shp[1][i-1] = 0.25*tn*(ss+tt+tt);
       shp[2][i-1] = 0.25*(1.+ss)*(1.+tt)*(-1.+ss+tt);
   }

   *xsj = (xs*yt-ys*xt)/4.0;

   xs = xs/ *xsj;
   xt = xt/ *xsj;
   ys = ys/ *xsj;
   yt = yt/ *xsj;
   *xsj = *xsj/4.0;

   for(i=5;i<=7;i=i+2) {
       ss = si[i-4-1]*s;
       tt = ti[i-4-1]*t;
       shp[0][i-1] =     -s*(1.+tt);
       shp[1][i-1] =    0.5*ti[i-4-1]*(1.-s*s);
       shp[2][i-1] =    0.5*(1.-s*s)*(1.+tt);
       shp[0][i] = -0.5*si[i-4-1]*(1.-t*t);
       shp[1][i] =   -t*(1.-ss);
       shp[2][i] =  0.5*(1.-ss)*(1.-t*t);
   }

   for(i=1; i<= 8; i++) {
      sn = yt*shp[0][i-1]-ys*shp[1][i-1];
      shp[1][i-1] = xs*shp[1][i-1]-xt*shp[0][i-1];
      shp[0][i-1] = sn;
   }

#ifdef DEBUG
       printf("\n");
       printf("In Qushp8 () : yt = %12.5f ys = %12.5f\n", yt, ys);
       printf("In Qushp8 () : xt = %12.5f xs = %12.5f\n", xt, xs);
       dMatrixPrint("shape func", shp, 3, 8);
#endif

   return(shp);
}

/*
 *  ===============================================
 *  Print DKT_PLATE Element Properties
 *  ===============================================
 */

#ifdef __STDC__
void print_property_plate(EFRAME *frp, int i)
#else
void print_property_plate(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[11].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Plate Thickness = %16.3e %s\n",
                           eap->work_section[11].value/eap->work_section[11].dimen->scale_factor,
                           eap->work_section[11].dimen->units_name);
	}
       break;
       case OFF:
        if( eap->work_material[0].value != 0.0 ) {
           printf("             ");
           printf("         : Young's Modulus =  E = %16.3e\n",
                            eap->work_material[0].value);
        }
        if( eap->work_material[2].value != 0.0 ) {
           printf("             ");
           printf("         : Yielding Stress = fy = %16.3e\n",
                            eap->work_material[2].value);
        }
        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[0].value != 0.0 ) {
           printf("             ");
           printf("         : Density         = %16.3e\n",
                            eap->work_material[5].value);
        }
        if( eap->work_section[11].value != 0.0 ) {
           printf("             ");
           printf("         : Plate Thickness = %16.3e\n",
                            eap->work_section[11].value);
        }
        break;
        default:
        break;
     }
}

ARRAY *sld08(p, isw)
ARRAY *p;
int isw;
{
    printf("ERROR >> ***In sld08() : elmt no =%3d : isw= %3d\n",p->elmt_no, isw);
    return(p);
}