/* *  =============================================================================  *  ALADDIN Version 1.0 : *       elmt_frame3d.c : Three Dimensional Frame Element *                                                                      *  Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin *  Institute for Systems Research,                                            *  University of Maryland, College Park, MD 20742                                    *                                                                      *  This software is provided "as is" without express or implied warranty. *  Permission is granted to use this software for any on any computer system *  and to redistribute it freely, subject to the following restrictions: *  *  1. The authors are not responsible for the consequences of use of *     this software, even if they arise from defects in the software. *  2. The origin of this software must not be misrepresented, either *     by explicit claim or by omission. *  3. Altered versions must be plainly marked as such, and must not *     be misrepresented as being the original software. *  4. This notice is to remain intact. *                                                                     *  -------------------------------------------------------------------  *  Three Dimensional Frame Element                                     *                                                                     *  Convention for Nodal Forces                                         *             +ve M     -  anticlockwise(RHT Rule)                     *             +ve X,Y,Z -  along +ve axis                              *  Convention for Member End  Forces                                   *             +ve M     -  Sagging Moments                            *             +ve SF    -  LHS upwards                                *             +ve AF    -  Tension(LHS outwards)                      *  -------------------------------------------------------------------  *                                                                      *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995 *  =============================================================================  */#include <math.h>#include "defs.h"#include "units.h"#include "matrix.h"#include "vector.h"#include "fe_database.h"#include "symbol.h"#include "fe_functions.h"#include "elmt.h"/*#define DEBUG*//* ============================================================== *//*   Element FRAME_3D                                             *//*   3D   Frame Element                                           *//*        Input Properties:                                       *//* ============================================================== *//*    p->work_material[0] = E;      p->work_material[1] = G;      p->work_material[2] = fy;      p->work_material[3] = ET;      p->work_material[4] = nu;      p->work_material[5] = density;      p->work_material[6] = fu;      p->work_section[0] = Ixx;      p->work_section[1] = Iyy;      p->work_section[2] = Izz;      p->work_section[3] = Ixy;      p->work_section[4] = Ixz;      p->work_section[5] = Iyz;      p->work_section[6] = weight;      p->work_section[7] = bf;      p->work_section[8] = tf;      p->work_section[9] = depth;                                        p->work_section[10] = area;      p->work_section[11] = plate_thickness;      p->work_section[12] = J;      p->work_section[13] = rT;      p->work_section[14] = width;      p->work_section[15] = tw;                                   *//* ============================================================== *//* macro for J constant */#define J_Const(x,y) ((1 - 0.63 * x/y)*( x * x* x* y/3))#ifdef __STDC__ARRAY *elmt_frame_3d(ARRAY *p, int isw)#elseARRAY *elmt_frame_3d(p, isw)ARRAY *p;int   isw;#endif{static double  nu;static QUANTITY  fy, G, E, ET, density;static double  Ixx, Iyy, Izz, Ixy, Ixz, Iyz, bf, tf, A, depth, weight, EA, EIzz, J, rT;double d1;double cs, sn,tn,xl,xx,yy,zz,vv,xn,xm,mbar;double **rot, **trot, **fr, **dlocal;int    i,j,k,l;DIMENSIONS  *dp_length, *dp_force, *dp_moment;DIMENSIONS  *dp_stress, *dp_degree, *dp_temperature;int          UNITS_SWITCH;#ifdef DEBUG       printf("*** Enter elmt_frame_3d() : isw = %4d\n", isw);#endif   H_Print = 0;   UNITS_SWITCH = CheckUnits();   switch(isw) {       case PROPTY:  /* MAT PROPS */          E.value       =  p->work_material[0].value;          fy.value      =  p->work_material[2].value;          ET.value      =  p->work_material[3].value;          nu            =  p->work_material[4].value;          density.value =  p->work_material[5].value;          if( UNITS_SWITCH == ON ) {             E.dimen       =  p->work_material[0].dimen;             fy.dimen      =  p->work_material[2].dimen;             ET.dimen      =  p->work_material[3].dimen;             density.dimen =  p->work_material[5].dimen;          }          /* (1)   check  poi_ratio value */          if( nu == 0.0 || nu > 0.5 ) {              printf("WARNING >> ... In 3d beam element() - frame_2d -  nu value = %9.4f,reset to 0.3 !\n", nu);              nu = 0.3;    /* default poi_ratio value */          }          /* (2)   calculate  G value */            /*        if(E.value/((1.0 - 2.0*nu)) != p->work_material[1].value) {              printf(" elmt_frame_3d(): WARNING: G is not equal to E/(1-2nu), check G for homogeneous material \n");              printf("                : ignore this message for non-homogeneous materials \n");          }*/          G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;          if(UNITS_SWITCH == ON)  G.dimen = E.dimen;          Ixx    = p->work_section[0].value;          Iyy    = p->work_section[1].value;          Izz    = p->work_section[2].value;          Ixy    = p->work_section[3].value;          Ixz    = p->work_section[4].value;          Iyz    = p->work_section[5].value;          weight = p->work_section[6].value;          bf     = p->work_section[7].value;          tf     = p->work_section[8].value;          depth  = p->work_section[9].value;          A      = p->work_section[10].value;          J      = p->work_section[12].value;          rT     = p->work_section[13].value;          EA     = E.value*A;          EIzz   = E.value*Izz;          /* (3) If J value not input, J calculated based on rectangular */          /*     section of size (bf x depth)                            */          if(J == 0.0 ) {             if(bf == 0.0 || depth == 0.0){                printf("WARNING >> Must give 'J' or ('width' & 'depth') to calculate stiffness");                exit(1);             }             /* Check bf < depth & cal J */             if(bf < depth )                  J = p->work_section[12].value = J_Const(bf,depth);              else                 J = p->work_section[12].value = J_Const(depth,bf);           }          break;       case CHERROR:            break;       case STRESS_LOAD:            break;       case STRESS_UPDATE:            break;       case PRESSLD:            break;       case STIFF:            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */            p->length.value = xl;            /* T matrix is made here: 12*12 size */            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);            rot = (double **) tmat(rot, 6, p);             if( UNITS_SWITCH==ON )   p->stiff->spColUnits[0].units_type = UNITS_TYPE;            p->stiff = beamst3d(p, p->stiff, EA, EIzz, E.value*Iyy, G.value*J ,xl, rot, p->size_of_stiff, p->dof_per_node);             MatrixFreeIndirectDouble(rot, p->size_of_stiff);            break;       case MASS_MATRIX:            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */            p->length.value = xl;            /* T matrix is made here: 12*12 size */            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);            rot = (double **) tmat(rot, 6, p);             /* Assemble Mass Matrix */             /* Calculate mbar =  mass/length                 */            /* in units of (kg/m) or (lbf*sec^2/in/in)       */            /* if no units, assume gravity g=9.80665 m/sec^2 */             if( weight != 0.0 )               mbar = weight/9.80665;             else               if( density.value > 0 )  mbar = A * density.value ;             else {                printf("\nError in input: Need density value to calculate mass matrix\n");                exit(1);             }	    /* Calculate radius of gyration , rT  --  m  ,  in   */            /* original version      :  rT = p->length.value/ 1.414; */            if( rT==0 && A!=0 )    rT = sqrt( J / A );            if( UNITS_SWITCH == ON )                 p->stiff->spColUnits[0].units_type = UNITS_TYPE;            p->stiff = beamms3d(p, p->stiff,p->type, mbar, xl, rT, rot, p->size_of_stiff, p->dof_per_node);            MatrixFreeIndirectDouble(rot, p->size_of_stiff);            break;       case STRESS:       case LOAD_MATRIX:            cs = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */            sn = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */            tn = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */            xl = sqrt(cs * cs + sn * sn + tn * tn); /* Calculate  Length */            p->length.value = xl;            /* T matrix is made here: 12*12 size */            rot = (double **) MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);            rot = (double **) tmat(rot, 6, p);             /* --------------------------- */            /* Output Stresses and Strains */            /* --------------------------- */            cs = rot[0][0];            rot[0][0] = 1.0;            if (UNITS_SWITCH==ON)  p->stiff->spColUnits[0].units_type = UNITS_TYPE;            p->stiff = beamst3d(p, p->stiff, EA, EIzz, E.value*Iyy,                       G.value*J, xl, rot, p->size_of_stiff, p->dof_per_node);             rot[0][0] = cs;            fr = MatrixAllocIndirectDouble(p->size_of_stiff,1); /* Here size_of_stiff = 12; i.e 12x1 mat */            for(l = 1; l<= p->nodes_per_elmt; l++) {                for(k = 1; k <= p->dof_per_node; k++) {                    j = (l-1)*p->dof_per_node + k;                    fr[j-1][0] = p->displ->uMatrix.daa[k-1][l-1];                }            }            dlocal = (double **)  dMatrixMult(rot, p->size_of_stiff, p->size_of_stiff, fr, p->size_of_stiff, 1);            fr = (double **)  dMatrixMultRep(fr, p->stiff->uMatrix.daa, p->size_of_stiff, p->size_of_stiff, dlocal, p->size_of_stiff, 1);            xx = 0.5 *(p->coord[0][0].value +p->coord[0][1].value);            yy = 0.5 *(p->coord[1][0].value +p->coord[1][1].value);            zz = 0.5 *(p->coord[2][0].value +p->coord[2][1].value);            if(H_Print == YES && isw == STRESS){               printf( "\n");               printf( "3D Frame Element   Stresses\n");               printf( "---------------------------\n");               H_Print = NO;            }            if(isw == STRESS && PRINT_STRESS == ON ) {               printf("\n");               printf("Elmt# %3d : ", p->elmt_no);               if(UNITS_SWITCH == ON)                  printf(": Coords (X,Y,Z)= (%8.3f %s,%8.3f %s,%8.3f %s)\n",                            xx/p->coord[0][0].dimen->scale_factor, p->coord[0][0].dimen->units_name,                            yy/p->coord[1][0].dimen->scale_factor, p->coord[1][0].dimen->units_name,                            zz/p->coord[2][0].dimen->scale_factor, p->coord[2][0].dimen->units_name);	       else                  printf(": Coords (X,Y,Z)= (%8.3f,%8.3f,%8.3f)\n", xx, yy, zz);               printf("\n");	    }            /*--------------------------------------------------------*/            /* nodal forces   & member end forces                     */            /*--------------------------------------------------------*/            if(p->elmt_load_ptr != NULL ) { /* calculate FEF */               printf("Fixed End Loads; \n");               p = sld05(p, STRESS);               /* Add FEF to joint p->[  ]  orces */                for(i = 1; i <= 12; i++)                  fr[i-1][0] = fr[i-1][0] - p->nodal_loads[i-1].value;            }            /* -------------------- */            /* Print Element forces */            /* -------------------- */            for(j = 1; j <= p->size_of_stiff; j++)                 p->nodal_loads[j-1].value = fr[j-1][0];            /* Assign element forces's units */            if( UNITS_SWITCH == ON ) {                if(UNITS_TYPE == SI) {                  dp_length = DefaultUnits("m");                  dp_force  = DefaultUnits("N");                }                if(UNITS_TYPE == US) {                  dp_length = DefaultUnits("in");                  dp_force  = DefaultUnits("lbf");                }                dp_moment = UnitsMult( dp_force, dp_length );                for(j=1;j<=3;j++) {		    UnitsCopy( p->nodal_loads[j-1].dimen, dp_force );		    UnitsCopy( p->nodal_loads[j-1+p->dof_per_node].dimen, dp_force );		    UnitsCopy( p->nodal_loads[j-1+3].dimen, dp_moment );		    UnitsCopy( p->nodal_loads[j-1+3+p->dof_per_node].dimen, dp_moment );		}            }            if(isw == STRESS && PRINT_STRESS == ON ) {	       switch(UNITS_SWITCH) {		 case ON:                  /* node_i */                  printf("            Fx1 = %13.5e %s\t Fy1 = %13.5e %s\t Fz1 = %13.5e %s\n",                         p->nodal_loads[0].value/dp_force->scale_factor, dp_force->units_name,                         p->nodal_loads[1].value/dp_force->scale_factor, dp_force->units_name,                         p->nodal_loads[2].value/dp_force->scale_factor, dp_force->units_name);                  printf("            Mx1 = %13.5e %s\t My1 = %13.5e %s\t Mz1 = %13.5e %s\n",                         p->nodal_loads[3].value/dp_moment->scale_factor, dp_moment->units_name,                         p->nodal_loads[4].value/dp_moment->scale_factor, dp_moment->units_name,                         p->nodal_loads[5].value/dp_moment->scale_factor, dp_moment->units_name);                  printf("\n");                  /* node_j */                  printf("            Fx2 = %13.5e %s\t Fy2 = %13.5e %s\t Fz2 = %13.5e %s\n",                        p->nodal_loads[6].value/dp_force->scale_factor, dp_force->units_name,                        p->nodal_loads[7].value/dp_force->scale_factor, dp_force->units_name,                        p->nodal_loads[8].value/dp_force->scale_factor, dp_force->units_name);                  printf("            Mx2 = %13.5e %s\t My2 = %13.5e %s\t Mz2 = %13.5e %s\n",                        p->nodal_loads[9].value/dp_moment->scale_factor, dp_moment->units_name,                        p->nodal_loads[10].value/dp_moment->scale_factor, dp_moment->units_name,                        p->nodal_loads[11].value/dp_moment->scale_factor, dp_moment->units_name);                  printf("\n");                  printf("            Axial Force : x-direction = %13.5e %s \n",                            -p->nodal_loads[0].value/dp_force->scale_factor,                             dp_force->units_name);                  printf("            Shear Force : y-direction = %13.5e %s \n",                             p->nodal_loads[1].value/dp_force->scale_factor,                             dp_force->units_name);                  printf("                        : z-direction = %13.5e %s \n",                             p->nodal_loads[2].value/dp_force->scale_factor,                             dp_force->units_name);                  printf("\n");                  break;                 case OFF:                  /* node_i */                  printf("            Fx1 = %13.5e\t Fy1 = %13.5e\t Fz1 = %13.5e\n",                           p->nodal_loads[0].value,                           p->nodal_loads[1].value,                           p->nodal_loads[2].value);                   printf("            Mx1 = %13.5e\t My1 = %13.5e\t Mz1 = %13.5e\n",                           p->nodal_loads[3].value,                           p->nodal_loads[4].value,                           p->nodal_loads[5].value);                   printf("\n");                   /* node_j */                  printf("            Fx2 = %13.5e\t Fy2 = %13.5e\t Fz2 = %13.5e\n",                           p->nodal_loads[6].value,                           p->nodal_loads[7].value,