elmt_frame3d.c

有限元程序

     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[1].dimen );
        UnitsSimplify( eap->work_section[2].dimen );
        UnitsSimplify( eap->work_section[12].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[1].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Inertia Iyy     = %16.3e %s\n",
                           eap->work_section[1].value/eap->work_section[1].dimen->scale_factor,
                           eap->work_section[1].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[12].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Torsional Constant J = %16.3e %s\n",
                           eap->work_section[12].value/eap->work_section[12].dimen->scale_factor,
                           eap->work_section[12].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[1].value != 0.0 ) {
           printf("             ");
           printf("         : Inertia Iyy     = %16.3e\n",
                            eap->work_section[1].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[12].value != 0.0 ) {
           printf("             ");
           printf("         : Torsional Constant J = %16.3e\n",
                            eap->work_section[12].value);
        }
        if( eap->work_section[10].value != 0.0 ) {
           printf("             ");
           printf("         : Area            = %16.3e\n",
                            eap->work_section[10].value);
        }
        break;
        default:
        break;
     }
}

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

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

/* printf(">>>>In sld05 : 3d elmt    \n");  */

       /* Initialize total load */

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

    switch(task){
        case PRESSLD: case STRESS:
             L = p->length.value;  elsptr  =  p->elmt_load_ptr;

             for(j=1; j<= elsptr->no_loads_faces; j++) {
                 elp  =  &elsptr->elib_ptr[j-1];

                 P =  elp->P.value;
                 a =  elp->a.value;
                 b =  elp->b.value;

                 /* Pt loads */

                 px =  elp->px.value;
                 py =  elp->py.value;
                 pz =  elp->pz.value;

                 /* moments */

                 mx =  elp->mx.value;
                 my =  elp->my.value;
                 mz =  elp->mz.value;

                 /* distributed loading */

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

                 if(a > L) /* error message */
                    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;
                    Z1 = bz * f2 * L;
                    MX1 = 0.0 ;
                    MY1 =-bz * f3 * L;
                    MZ1 = by * f3 * L;

                    X2 = bx * f4 * L;
                    Y2 = by * f5 * L;
                    Z2 = bz * f5 * L;
                    MX2 = 0.0 ;
                    MY2 =-bz * f6 * L;
                    MZ2 = by * f6 * L;

                    if (task == STRESS){
                        /* +ve  simply support moment at center */
                        if(b==L && a== 0.0) {/* udl acting on entire length */
                           MCZ = -by * (L * L)/8;   
                           MCY = -bz * (L * L)/8;   
                        }
                    else { /* approximate mom at center */
                        MCZ = -by *(b-a)* (L - (a+b)/2)/2;   
                        MCY = -bz *(b-a)* (L - (a+b)/2)/2;   
                    }
                 }
               }      /* end of dist loading */

               else {                 /* Concentrated Loading Condition */
                   /* shape functions */
                  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;

                  if(my != 0.0 || mz != 0.0){
                     /* derivatives of shape function */  
                     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;
                  Z1 = pz * f2 - my *  df2x;
                  MX1 =       mx * f1 ;
                  MY1 =-pz * f3 + my *  df3x;
                  MZ1 = py * f3 + mz *  df3x;

                  X2 = px * f4 + 0;
                  Y2 = py * f5 + mz *  df5x;
                  Z2 = pz * f5 - my *  df5x;
                  MX2 =       mx * f4 ;
                  MY2 =-pz * f6 + my *  df6x;
                  MZ2 = py * f6 + mz *  df6x;

                  if(task == STRESS) { /* +ve  simply support moment at center */
                     MCY = -pz * (L -a)/2 ;   
                     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] + Z1;
                 load[3] = load[3] + MX1;
                 load[4] = load[4] + MY1;
                 load[5] = load[5] + MZ1;
                 load[6] = load[6] + X2;
                 load[7] = load[7] + Y2;
                 load[8] = load[8] + Z2;
                 load[9] = load[9] + MX2;
                 load[10] = load[10] + MY2;
                 load[11] = load[11] + MZ2;

                 if(task == STRESS) { /* mid pt moment */
                    MCYT = MCYT+ MCY;
                    MCZT = MCZT+ MCZ;
                 }
              }
              break;
         default:
              break;
    }
  
    /* Rotate to Global */
 
    rmat = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
    rmat = (double **) tmat(rmat,6,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 = p->nodal_loads[inc-1].value + rt[inc-1][j-1]* (double) load[j-1];
    }

    MatrixFreeIndirectDouble(rmat, p->size_of_stiff);
    MatrixFreeIndirectDouble(rt, p->size_of_stiff);

    return(p);
}