elmt_plate.c

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

                   kk = p->dof_per_node*(i-1) + 1;
                   for(j = 1; j <= p->dof_per_node; j++) {
                       k  = p->dof_per_node*(i-1) + j;
                       if( k <= kk) {
                           UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) );
                           UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[0]) );
                       }
                       if(k > kk) {
                           UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[1]) );
                           UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[1]) );
                       }
                   }
               }

               free((char *) dimen1->units_name);
               free((char *) dimen1);
               free((char *) dimen2->units_name);
               free((char *) dimen2);
            }
     
            break;
       case STRESS:   /* Compute And Output The Element Variables */

            jacqud(p->coord,aa,bb,cc,dd,ee);
            shp = qushp8(0.,0.,shp,p->coord,&xsj);
            bm  = dktqbm(shp,bm,aa,bb,cc,dd,ee);

            for(i=1; i<=3; i++) {
                eps[i-1] = 0.0;
                for(j=1; j <= 4; j++) {
                    for(k=1; k <= 3; k++) 
			eps[i-1] += bm[i-1][(int) (3*(j-1)+k-1)]*p->displ->uMatrix.daa[k-1][j-1];
                }
            }

            sig[0] = (p->work_section[12].value*eps[0])+(p->work_section[13].value*eps[1]);
            sig[1] = (p->work_section[13].value*eps[0])+(p->work_section[12].value*eps[1]);
            sig[2] =  p->work_section[14].value*eps[2]/2.0;

            xx = 0.25*(p->coord[0][0].value+p->coord[0][1].value+p->coord[0][2].value+p->coord[0][3].value);
            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(UNITS_TYPE) {
	            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) */
/* ============= */

#ifdef OLD_VERSION
/***********************/
int jacqud(x,ndm,aa,bb,cc,dd,ee)
double **x;
double aa[5],bb[5],cc[5],dd[5],ee[5];
int ndm;
{
int i,k;
double b,c,sql;

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

       b = x[3][k]-x[3][i]; 
       c = x[1][i]-x[1][k];
       sql = b*b+c*c;

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

   return;
}
#endif
/***********************/

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

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