elmt_frame2d.c

有限元程序

/*
 *  ===============================================================
 *  rotate () : Rotate element to global coordinate frame
 *
 *  Input  : double **s   -- pointer to stiffness/mass matrix.
 *         : double cs    -- cosine of the element orientation.
 *         : double sn    -- sine of the element orientation.
 *         : int nst      -- 
 *         : int ndf      -- maximum d.o.f. at any node.
 *  Output : MATRIX *     -- pointer to rotated matrix.
 *  ===============================================================
 */

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

   /* [a] : Check to see if element orientation is horizontal */

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

   /* [b] : Rotate matrix to global coordinate frame */

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


/*
 *  ===============================================================
 *  print_property_frame_2d() : print FRAME_2D Element Properties
 *
 *  Input  : EFRAME *frp  --
 *         : int i        -- 
 *  Output : void 
 *  ===============================================================
 */

#ifdef __STDC__
void print_property_frame_2d(EFRAME *frp, int i)
#else
void print_property_frame_2d(frp, i)
EFRAME    *frp;
int          i;                 /* elmt_attr_no */
#endif
{
int     UNITS_SWITCH;
ELEMENT_ATTR    *eap;

#ifdef DEBUG
       printf("*** Enter print_property_frame_2d()\n");
#endif

     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[2].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Inertia Izz     = %16.3e %s\n",
                           eap->work_section[2].value/eap->work_section[2].dimen->scale_factor,
                           eap->work_section[2].dimen->units_name);
	}
	if( eap->work_section[10].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Area            = %16.3e %s\n",
                           eap->work_section[10].value/eap->work_section[10].dimen->scale_factor,
                           eap->work_section[10].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[2].value != 0.0 ) {
           printf("             ");
           printf("         : Inertia Izz     = %16.3e\n",
                            eap->work_section[2].value);
        }
        if( eap->work_section[10].value != 0.0 ) {
           printf("             ");
           printf("         : Area            = %16.3e\n",
                            eap->work_section[10].value);
        }
        break;
        default:
        break;
     }
#ifdef DEBUG
       printf("*** Leave print_property_frame_2d()\n");
#endif
}


/* 
 *  ====================================================================
 *  sld07() : Compute equivalent nodal loads for 2d frame element
 * 
 *  Settings for the task flag are:
 * 
 *      task         Action
 *      ======================================
 *      PRESSLD      Pressure loading.
 *      STRESS       Streess  loading.
 * 
 *  Input  : ARRAY *  -- pointer to working array data structure.
 *         : int task -- flag for task to be performed.
 *  Output : ARRAY *  -- pointer to working array data structure.
 *  ====================================================================
 */ 

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

    /* [a] : 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()\n" );  
                printf("ERROR >> 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;
       default:
            break;
    }

    /* [b] : 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];
         }
    }

    /* [c] : Mid point moment */
 
    MatrixFreeIndirectDouble(rmat, p->size_of_stiff);
    MatrixFreeIndirectDouble(rt, p->size_of_stiff);

    return(p);
}