elmt_psps.c

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

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

}