elmt_psps.c

有限元分析源代码

    t[0] =  0.5; t[1] =  0.5;
    t[2] = -0.5; t[3] = -0.5;

  switch(nodes_per_elmt) { 
    case 3: 
    case 4:

    /* form 4-node quadrilateral shape function                  */
    /* shape function:                  Ni = shape[2][i]         */
    /*                                  node no. i = 1, .. 4     */
    /* derivatives of shape functions:  dNi/d(ss) = shape[0][i]  */
    /*                                  dNi/d(tt) = shape[1][i]  */

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

    /* form triangle by adding third and fourth node together */

    if(nodes_per_elmt == 3) { 
        shp[2][2] = shp[2][2] + shp[2][3];
        for(i = 0; i <= 1; i++)
            shp_temp[i][2] = shp_temp[i][2] + shp_temp[i][3];
    }

     /* construct jacobian matrix and its determinant */

     for(i = 1; i <= no_dimen; i++)      /* no_dimen = 2 */
         for(j = 1; j <= no_dimen; j++) {
             xs[i-1][j-1] = 0.0;
             for(k = 1; k <= nodes_per_elmt; k++)
                 xs[i-1][j-1] = xs[i-1][j-1] + coord[i-1][k-1].value*shp_temp[j-1][k-1];
         }

     *jacobian = xs[0][0] * xs[1][1] - xs[0][1] *xs[1][0];

    if(flg == TRUE) return;

    /* compute Jacobain inverse matrix */

    sx[0][0] = xs[1][1]/ *jacobian;
    sx[1][1] = xs[0][0]/ *jacobian;
    sx[0][1] = - xs[0][1]/ *jacobian;
    sx[1][0] = - xs[1][0]/ *jacobian;


    /* form global derivatives */

    /* save dNi/dx, dNi/dy into shp[2][node] */
  
    for(i = 1; i <= nodes_per_elmt; i++){
        shp[0][i-1] = shp_temp[0][i-1]*sx[0][0] + shp_temp[1][i-1]*sx[0][1];
        shp[1][i-1] = shp_temp[0][i-1]*sx[1][0] + shp_temp[1][i-1]*sx[1][1];
    
    }
    break;

    case 5: case 6: case 7: case 8:
       shp0(ss, tt, shp, node_connect, nodes_per_elmt);
    break;

    default:
    break;
  }
#ifdef DEBUG

   printf(" in shape () \n");
    for(j = 1; j <= nodes_per_elmt; j++){
            printf("   dN/dx[%d] = %lf ", j , shp[0][j-1]);
            printf("   dN/dy[%d] = %lf ", j , shp[1][j-1]);
            printf("       N[%d] = %lf \n", j , shp[2][j-1]);
   }
   printf(" leaving shape () \n");
#endif

   MatrixFreeIndirectDouble(shp_temp, 2);
    return;
}

/* ================================ */
/* shap0                            */
/* ================================ */

int shp0(s,t,shp,ix,nel)
double s,t,shp[4][10];
int  *ix;
int nel;
{
double s2, t2;
int i,j,k,l;

    s2 = (1- s * s)/2;
    t2 = (1 - t * t)/2;

    for(i = 5; i<= 9 ;i++)
    for(j =1;j<= 3 ; j++)
        shp[j][i] = 0.0;

    /* midside nodes (serendipty)  */

    if(ix[5] == 0) goto NEXT1;
        shp[1][5] = -s*(1-t);
        shp[2][5] = -s2;
        shp[3][5] =  s2*(1-t);

    NEXT1:
        if(nel < 6) goto NEXT7;
        if(ix[6] == 0) goto NEXT2;

        shp[1][6] = t2;
        shp[2][6] = - t *(1+s);
        shp[3][6] =   t2*(1+s);

   NEXT2:
        if(nel < 7) goto NEXT7;
        if(ix[7] == 0) goto NEXT3;

        shp[1][7] = -s*(1 + t);
        shp[2][7] = s2;
        shp[3][7] = s2*(1+t);

   NEXT3:
        if(nel < 8) goto NEXT7;
        if(ix[8] == 0) goto NEXT4;

        shp[1][8] = -t2;
        shp[2][8] = - t *(1-s);
        shp[3][8] = t2*(1-s);

     /* interior node (langragian) */
   NEXT4:
        if(nel < 9) goto NEXT7;
        if(ix[9] == 0) goto NEXT7;

        shp[1][9] = -4 * s *  t2;
        shp[2][9] = -4 * t * s2;
        shp[3][9] =  4 * s2 * t2;

     /* correct edge nodes for interior node(langrangian) */

        for(j=1;j<=3;j++) {
            for(i=1;i<=4;i++)
                shp[j][i] = shp[j][i] - 0.25 * shp[j][9];
                for(i=5;i<=8;i++)
                if(ix[i] != 0)
	           shp[j][i] = shp[j][i] - .5 *shp[j][9];
        }

     /* correct corner nodes for presence of mideside nodes */
   NEXT7:
        k = 8;
        for(i=1;i<=4;i++){
            l = i +4;
            for(j=1;j<=3;j++)
                shp[j][i] = shp[j][i] - 0.50 *(shp[j][k] + shp[j][l]);
                k = 1;
        }

        return(1);

}


/* ======================== */
/* gauss  integration       */
/* ======================== */

#ifdef __STDC__
gauss(double *sg, double *ag, int lt)
#else
gauss(sg, ag, lt)
double *sg, *ag;
int lt;
#endif
{
double t;
int    i;

#ifdef DEBUG
   printf(" ******enter gauss() \n");
   printf("       no_integ_pts = %d \n", lt);
#endif

switch(lt) {
  case 1:
         sg[1] = 0.0;
         ag[1] = 2.0;
         break;
  case 2:
         sg[1] = -1/sqrt(3.);
         sg[2] =  -sg[1];
         ag[1] = 1.0; 
         ag[2] = 1.0; 
         break;
  case 3:
         sg[1] = -sqrt(0.6);
         sg[2] = 0.0;
         sg[3] = -sg[1];

         ag[1] = 5.0/9.0; 
         ag[2] = 8.0/9.0; 
         ag[3] = ag[1]; 
         break;
  case 4:
         t = sqrt(4.8);
         sg[1] =  sqrt((3+t)/7);
         sg[2] =  sqrt((3-t)/7);
         sg[3] = -sg[2];
         sg[4] = -sg[1];
         t  = 1/3/t;
         ag[1] = 0.5-t; 
         ag[2] = 0.5+t; 
         ag[3] = ag[2]; 
         ag[4] = ag[1]; 
         break;
  default:
         break;
  }

#ifdef DEBUG
      printf(" In gauss(): \n");
   for (i = 1; i <= lt; i++) {
      printf("           : integ_coord[%d] = %lf \n",i, sg[i]);
      printf("           : weight[%d]      = %lf \n",i, ag[i]);
   }
   printf(" ******leaving gauss() \n");
#endif

  return (1);
}

/*====================================================*/
/* pgauss(no_integ_pts,lint,sg,tg,wg)                 */
/* input:                                             */
/*        no_integ_pts   no of gaussian integ pts     */
/*                       in one-direction             */
/* output:                                            */
/*        lint           no_integ_pts*no_integ_pts    */
/*        sg             coordinates in xi direction  */
/*        tg             coordinates in eta direction */
/*        wg             weight coeffcient            */
/*====================================================*/

int pgauss(l, lint, r, z, w)
int l, *lint;
double r[], z[], w[];
{
int i, j, k;
double g, h;
double g4[5], h4[5], lr[10], lz[10], lw[10];

lr[1] = lr[4] =lr[8] = -1;
lr[2] = lr[3] =lr[6] =  1;
lr[5] = lr[7] =lr[9] =  0;

lz[1] = lz[2] =lz[5] = -1;
lz[3] = lz[4] =lz[7] =  1;
lz[6] = lz[8] =lz[9] =  0;

lw[3] = lw[4] = lw[1] = lw[2] = 25;
lw[5] = lw[6] = lw[7] = lw[8] = 40;
lw[9] = 64;

    if(l < 0) {
       *lint = 3;
       g = sqrt((double) 1.0/3.0);
       h = sqrt((double) 2.0/3.0);
       r[1] = -h; r[2] =  h; r[3] = 0;
       z[1] = -g; z[2] = -g; z[3] = g;
       w[1] =  1; w[2] =  1; w[3] = 2;

       return (1);
    }

    *lint = l*l;
    switch (l) {
	case 1: /*  1 x 1 integration */
	     r[0] = 0.0;
	     z[0] = 0.0;
	     w[0] = 4.0;
	     break;
        case 2: /*  2 x 2 integration */
             g = 1.0/sqrt((double) 3.0);
             for(i = 1; i<= 4; i++) {
                 r[i-1] = g * lr[i];
                 z[i-1] = g * lz[i];
                 w[i-1] = 1.0;
             }
             break;
        case 3: /* 3 x 3 integration */
             g = sqrt((double) 0.6);
             h = 1.0/81.0;
             for(i = 1; i<= 9; i++) {
                 r[i-1] = g * lr[i];
                 z[i-1] = g * lz[i];
                 w[i-1] = h * lw[i];
             }
             break;
        case 4: /* 4 x 4 integration */
             g = sqrt((double) 4.8);
             h = sqrt((double) 30.0)/36;
             g4[1] = sqrt((double) (3+g)/7.0);
             g4[4] = -g4[1];
             g4[2] = sqrt((double) (3-g)/7.0);
             g4[3] = -g4[2];
             h4[1] = 0.5 - h;
             h4[2] = 0.5 + h;
             h4[3] = 0.5 + h;
             h4[4] = 0.5 - h;
             i = 0;
             for(j = 1; j<= 4; j++) {
                 for(k = 1; k<= 4; k++) {
                     i = i +1;
                     r[i-1] = g4[k];
                     z[i-1] = g4[j];
                     w[i-1] = h4[j]* h4[k];
                 }
             }
             break;
    }

    return(1);
}


/* ========================================== */
/* pstres                                     */
/* ========================================== */

int pstres(sig,sg4,sg5,sg6)
double *sig;
double sg4,sg5,sg6;
{
    printf("******In function pstres\n ");
    return(1);
}

/*=============================================*/
/* material matrix                             */
/*=============================================*/
double **MATER_MAT_PLANE(E, nu)
double E, nu;
{
 
double **m1;
double temp;

    m1 = MatrixAllocIndirectDouble(3, 3);
    
    temp =  E/(1.0-nu*nu);

    m1[0][0]= m1[1][1] = temp;
    m1[0][1]= m1[1][0] = nu*temp; 
    m1[2][2]= E/2.0/(1.0+nu); 

    m1[0][2] = m1[2][0] = m1[1][2] = m1[2][1] = 0.0;

    return(m1);
}


/* =================================================================== */
/* Calculate Equivalent Nodal Loads                                    */
/*                                                                     */
/* These element level procedures determine Eqvivalent Nodal Loads due */
/* to element loadings. The loads are added directly to nodal loads.   */
/*                                                                     */
/*  Negative Eqv (= Fixed end forces) added to local elmt forces       */
/*   Sign Convention:  P, by,bx,etc positive along Y,X axis            */
/*   sld01 Equivalent Loading Procedure for 2 D finite element         */
/* =================================================================== */

ARRAY *sld01( p,task)
ARRAY *p;
int task;
{
ELEMENT         *el;  
ELEMENT_LOADS  *ell; 
ELOAD_LIB      *elp; 

double shp[4][10],sg[17],tg[17],wg[17],pq[4],dj[4],xyp[4],xye[4][5];
double h, dv,dv1,dv2,ax,ay,jacobain;
int mash,gash,gish,gosh,nash,l,lint,nsfr,lne,ncord,nln;
int i,j,lnew, k,igauss;
int no_integ_pts;

    switch(task){
        case PRESSLD:
             nsfr  = 2;/* conter for variation type; para =2 */
             lne   = 4;/* no of element face */
             ncord = 3;
             nln = nsfr + 1;
             ell  =  p->elmt_load_ptr;
             for(lnew=1;lnew<=lne; lnew++){
                 h = -1;
                 for (j=1; j<=nln; j++){
                      elp  =  &ell->elib_ptr[lnew-1];
                      mash =   elp->nopl[j];
                      for(k=1; k<=ncord; k++)
                          xye[k][j] = p->coord[mash][k].value;
                 }

            if(p->nodes_per_elmt == 4) no_integ_pts  =  2;    /* 4-node element */
            if(p->nodes_per_elmt == 8) no_integ_pts  =  3;    /* 8-node element */
            l = no_integ_pts;

                 pgauss(l,lint,sg,tg,wg);

                 for(igauss=1;igauss<=lint;igauss++){
                     shape(sg[igauss],h,p->coord,shp,jacobain,p->no_dimen,p->nodes_per_elmt, p->node_connect,0);
                     for(i = 1; i <= ncord; i++){
                         gash = 0;
                         gosh = 0;
                         gish = 0;
                         for (k=1; k<=nln; k++){
                              if(i<=2){;
                                 gosh = gosh +   elp->pr->uMatrix.daa[k][i] * shp[3][k];
                                 gish = gish + xye[i][k] * shp[1][k];
                              }
                         }
                         xyp[i] = gash;
                         if(i<=2){
                            pq[i] = gosh;
                            dj[i] = gish;
                         }

                     }
                     dv = wg[igauss];
                     if(ncord == 3) dv = dv + xyp[3];
                     if(p->nodes_per_elmt == AXISYM) dv = 6.283185 * xyp[1];
                        dv1 = dv + dj[1];
                        dv2 = dv + dj[2];
                        ay = dv1*pq[1] + dv2*pq[2];
                        ax = dv1 * pq[2] - dv2 * pq[1];
                        for (i=1; i<=nln; i++){
                             nash = p->dof_per_node * ( ell->elib_ptr[lnew -1].nopl[i] - 1)  + 2; 
                             mash = nash -1;
                             F[mash] = F[mash] + shp[3][i] * ax;
                             F[nash] = F[nash] + shp[3][i] * ay;
                        }
                 }
            }
            break;
       default:
            break;
    }

    return(p);

}