elmt_frame2d.c

有限元分析源代码

	     break;
	default:
             FatalError("In elmt_frame2d() : beamms() : Type of Mass Matrix Undefined",(char*)NULL);
	     break;
    }

    /* ==================== */
    /* Initial units buffer */
    /* ==================== */

    if( CheckUnits() == ON ) {
       if(UNITS_TYPE == SI) {
          d1 = DefaultUnits("Pa");
          d2 = DefaultUnits("m");
       }
       else {
          d1 = DefaultUnits("psi");
          d2 = DefaultUnits("in");
       }
       d3 = DefaultUnits("sec");

       UnitsMultRep( &(s->spColUnits[0]), d1, d2 );
       UnitsCopy( &(s->spColUnits[1]), &(s->spColUnits[0]) );
       UnitsMultRep( &(s->spColUnits[2]), &(s->spColUnits[0]), d2 );

       UnitsPowerRep( &(s->spRowUnits[0]), d3, 2.0, NO );
       UnitsCopy( &(s->spRowUnits[1]), &(s->spRowUnits[0]) );
       UnitsMultRep( &(s->spRowUnits[2]), d2, &(s->spRowUnits[0]) );

       UnitsCopy( &(s->spColUnits[3]), &(s->spColUnits[0]) );
       UnitsCopy( &(s->spColUnits[4]), &(s->spColUnits[1]) );
       UnitsCopy( &(s->spColUnits[5]), &(s->spColUnits[2]) );

       UnitsCopy( &(s->spRowUnits[3]), &(s->spRowUnits[0]) );
       UnitsCopy( &(s->spRowUnits[4]), &(s->spRowUnits[1]) );
       UnitsCopy( &(s->spRowUnits[5]), &(s->spRowUnits[2]) );

       free((char *) d1->units_name);
       free((char *) d1);
       free((char *) d2->units_name);
       free((char *) d2);
       free((char *) d3->units_name);
       free((char *) d3);
    }

    /* ==================== */
    /* Rotate mass matrix   */
    /* ==================== */

    s->uMatrix.daa = (double **) rotate(s->uMatrix.daa,cs,sn,nst,ndf);

#ifdef DEBUG
       MatrixPrint(s);
       printf("*** leaving beamms()  \n");
#endif

    return(s);
}

/* ====================== */
/* Rotate                 */
/* ====================== */

#ifdef __STDC__
double **rotate(double **s, double cs, double sn, int nst, int ndf)
#else
double **rotate(s, cs, sn, nst, ndf)
double **s;
int  nst, ndf;
double cs, sn;
#endif
{
int    i,j,n;
double t;

   if(cs == 1.0)
      return(s);

   for(i = 0; i < nst; i = i + ndf){
       j = i+1;
       for(n = 0; n < nst;n++){
           t = s[n][i] * cs - s[n][j] * sn;
           s[n][j] = s[n][i] * sn + s[n][j] *cs;
           s[n][i] = t;
       }
   }

   for(i = 0; i < nst; i=i+ndf){
       j = i+1;
       for(n = 0; n < nst; n++){
           t = cs*s[i][n] - sn*s[j][n];
           s[j][n] = sn * s[i][n] + cs * s[j][n];
           s[i][n] = t;
       }
   }
 
   return(s);
}


/* =================================================== */
/* Equivalent Loading Procedure for 2 D frame element  */
/* =================================================== */

#ifdef __STDC__
ARRAY *sld07(ARRAY *p, int task)
#else
ARRAY *sld07(p,task)
ARRAY *p;
int task;
#endif
{
ELEMENT_LOADS *elsptr;
ELOAD_LIB        *elp;
double  P, a ,b ;
double  L, load[8];
double  px,py,pz, mx,my,mz, bx,by,bz, ze,ze2,ze3;
double  X1,Y1,MZ1,
        X2,Y2,MZ2, 
        MCZ, MCZT;
double  f1,f2,f3,f4,f5,f6, Z1, Z2,
        df1x, df3x, df5x,df2x, df6x;
double  **rmat,**rt;
int     *da;
int     inc,j, i;

    /* Initialize total load */

    for(inc=1; inc<=7;inc++)
        load[inc-1] = 0.0;
    MCZT = 0.0;

    switch(task){
        case PRESSLD:
        case STRESS:
             L       = (double) p->length.value;  
             elsptr  =  p->elmt_load_ptr;
             for (j=1; j<= elsptr->no_loads_faces;j++) {
                  elp = &elsptr->elib_ptr[j-1];

             /* Pt loads */

             P = elp->P.value;
	     a = elp->a.value;
	     b = elp->b.value;
             px = elp->px.value;
	     py = elp->py.value;
       
             /* moments */

             mz = elp->mz.value;

             /* distributed loading */

             bx =  elp->bx.value;
             by =  elp->by.value;

             if(a > L)
                printf(">>ERROR in sld; Elmt Load dist. 'a' > Elmt Length; El_no= %d\n",p->elmt_no);  

             if (elp->type == -1) { /* Distributed loading Condition */
                inc = 0;
                /* set default values */
                if(b == 0.0) b = L; /* dist loading acts on entire length */

                /* first calc f(b) */

                ze = b/L;    
SHP_START:
                ze2 = ze * ze; ze3 = ze2 * ze;
                f1 = 1   -  ze2/2;
                f2 = ze3*ze/2  - ze3 + ze ;
                f3 = (ze3*ze/4 - 2*ze3/3 + ze2/2) * L;
                f4 = ze2/2 ;
                f5 = -ze *ze3/2 +  ze3;
                f6 = (ze3*ze/4  - ze3/3) * L;
                inc++;
                 
                if( inc == 1) {
                   /* temp hold of values f(b) */
                   X1 = f1; Y1 = f2;  Z1 = f3; 
                   X2 = f4; Y2 = f5;  Z2 = f6; 
                   ze = a/L;
                   goto SHP_START;
                }
                else{
                   /* f() = f(b) - f(a)  */
                   f1 = X1 - f1; f2 =  Y1 - f2;f3 = Z1 - f3; 
                   f4 = X2 - f4; f5 =  Y2 - f5;f6 = Z2 - f6; 
                }

                X1 = bx * f1 * L;
                Y1 = by * f2  * L;
                MZ1 = by * f3 * L;

                X2  = bx * f4 * L;
                Y2  = by * f5 * L;
                MZ2 = by * f6 * L;

                /* +ve  simply support moment at center */

                if (task == STRESS){

                   if(b==L && a== 0.0)  /* udl acting on entire length */
                      MCZ = -by * (L * L)/8;   
                   else   /* approximate mom at center */
                    MCZ = -by *(b-a)* (L - (a+b)/2)/2;   
                }

            } /* end of dist loading */

	    /* Concentrated Loading Condition */
            /* distributed body force         */

            else {
              /* shape functions */
              /* =========================================*/
              /* according to :                           */
              /* "Structural Dynamics by Finite Elements" */
              /* W. Weaver, Jr. and P. R. Johnson         */
              /* Chapter 6: Space Frame                   */
              /* =========================================*/


              ze = a/L;      ze2 = ze * ze;        ze3 = ze2 * ze;

              f1 =     1  -  ze;
              f2 =     2 * ze3  -  3 * ze2 + 1;
              f3 =    (ze3 - 2 * ze2 + ze) * L;
              f4 =     ze ;
              f5 =    -2 * ze3  +  3 * ze2;
              f6 =    (ze3 - ze2 ) * L;

              /* derivatives of shape function */  

              if(mz != 0.0){
                 df2x =  6 *( ze2  -  ze) / L;
                 df3x =  3 *ze2  - 4*ze  + 1;
                 df5x = - df2x; 
                 df6x =  3 *ze2  - 2*ze;
              }

              X1 = px * f1 + 0;
              Y1 = py * f2 + mz *  df2x;
              MZ1 = py * f3 + mz *  df3x;

              X2 = px * f4 + 0;
              Y2 = py * f5 + mz *  df5x;
              MZ2 = py * f6 + mz *  df6x;

              /* +ve  simply support moment at center */
              if(task == STRESS) 
                 MCZ = -py * (L -a)/2 ;   
                  
            }

            /* Add Contributation to Total Equivalent Load */

               load[0] = load[0] + X1;
               load[1] = load[1] + Y1;
               load[2] = load[2] + MZ1;
               load[3] = load[3] + X2;
               load[4] = load[4] + Y2;
               load[5] = load[5] + MZ2;
	       
               if(task == STRESS)
	          MCZT = MCZT+ MCZ;
            }
            break;
       case STRESS_LOAD:
            for(i = 1; i<= 6; i++)
            load[i-1] = p->nodal_loads[i-1].value ;
            break;
       default:
            break;
    }

    /* Rotate Local forces to Global forces */
 
    rmat = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    rmat = (double **) tmat(rmat,4,p);
    rt   = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    for( i=1 ; i<=p->size_of_stiff ; i++ )
       for( j=1 ; j<=p->size_of_stiff ; j++ )
          rt[j-1][i-1] = rmat[i-1][j-1];

    for (inc=1; inc<=p->size_of_stiff; inc++){
         p->nodal_loads[inc-1].value = 0.0;
         for (j=1; j<=p->size_of_stiff; j++) {
              p->nodal_loads[inc-1].value += rt[inc-1][j-1]* (double) load[j-1];
         }
    }
    /* mid pt moment */
/*
    if (task == STRESS) {
        p->nodal_loads[7].value = MCZT; 
    }
*/
 
    MatrixFreeIndirectDouble(rmat, p->size_of_stiff);
    MatrixFreeIndirectDouble(rt, p->size_of_stiff);

    return(p);
}