elmt_frame3d.c

有限元程序

            /* Assign element forces's units */

            if( UNITS_SWITCH == ON ) {
                if( CheckUnitsType() == SI) {
                    dp_length = DefaultUnits("m");
                    dp_force  = DefaultUnits("N");
                }
                if( CheckUnitsType() == 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 ) {
	       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,
                           p->nodal_loads[8].value); 
                  printf("            Mx2 = %13.5e\t My2 = %13.5e\t Mz2 = %13.5e\n",
                           p->nodal_loads[9].value,
                           p->nodal_loads[10].value,
                           p->nodal_loads[11].value); 
                  printf("\n"); 

                  printf("            Axial Force : x-direction = %13.5e \n", -p->nodal_loads[0].value);
                  printf("            Shear Force : y-direction = %13.5e \n",  p->nodal_loads[1].value);
                  printf("                        : z-direction = %13.5e \n",  p->nodal_loads[2].value);
                  printf("\n"); 
                 break;
                 default:
                 break;
               }
	    }

            if(isw == LOAD_MATRIX ) {
                for(j = 1; j <= p->size_of_stiff; j++) 
                    fr[j-1][0] = 0.0;
                for( i=1 ; i<=p->size_of_stiff; i++ ) {
                    for( j=1 ; j<=p->size_of_stiff; j++ ) {
                        fr[i-1][0] += rot[j-1][i-1]*p->nodal_loads[j-1].value;
                    }
                }
                for(j = 1; j <= p->size_of_stiff; j++) 
                    p->nodal_loads[j-1].value = fr[j-1][0];
            }

            if( UNITS_SWITCH==ON ) {
                free((char *) dp_length->units_name);
                free((char *) dp_length);
                free((char *) dp_force->units_name);
                free((char *) dp_force);
                free((char *) dp_moment->units_name);
                free((char *) dp_moment);
            }

            MatrixFreeIndirectDouble(fr, p->size_of_stiff);
            MatrixFreeIndirectDouble(dlocal, p->size_of_stiff);
            MatrixFreeIndirectDouble(rot, p->size_of_stiff);
            break;
      default:
            break;
   }

   return(p);
}


/* =============================================================== */
/* 3D Frame Stiffness: This function returns a 12x12 rotated       */
/* stiffness matrix for 3D frame elmts.                            */
/*                                                                 */
/* Input:                                                          */
/*   ea->EA ; s->matrix pointer; eiz->EI_zaxis;  eiy->EI_yaxis;    */
/*   gj->G.J; xl->length       ; rot->rotation_matrix(12x12)       */
/*   nst->num_of_deg_of free per element; ndf->dof's per node;     */
/* =============================================================== */

#ifdef __STDC__
MATRIX *beamst3d(ARRAY *p, MATRIX *s, double ea, double eiz, double eiy, double gj, double xl, double **rot, int nst, int no_dof)
#else
MATRIX *beamst3d(p,s, ea, eiz, eiy, gj, xl, rot, nst, no_dof)
ARRAY  *p;
MATRIX *s;
double **rot;
double ea, eiz,eiy,gj, xl;
int    nst, no_dof;
#endif
{
int    i, j, k,ii,jj,kk;
double t;
DIMENSIONS     *d1, *d2;

    i  = no_dof + 1;
    j  = no_dof + 2; 
    k  = no_dof + 3;
    ii = no_dof + 4;
    jj = no_dof + 5;
    kk = no_dof + 6;

    t = ea/xl; 
#ifdef DEBUG

    printf(" in beamst3d() : \n");
    printf("               : EA = %lf, xl = %lf \n", ea, xl);
    printf("               : t = %lf\n", t);
#endif

    s->uMatrix.daa[0][0] =  t;
    s->uMatrix.daa[i-1][i-1] =  t;
    s->uMatrix.daa[0][i-1] = -t;
    s->uMatrix.daa[i-1][0] = -t;

    t = 12 * eiz/xl/xl/xl ;

    s->uMatrix.daa[1][1] = t;
    s->uMatrix.daa[j-1][j-1] = t;
    s->uMatrix.daa[1][j-1] = -t;
    s->uMatrix.daa[j-1][1] = -t;

    t = 12 * eiy/xl/xl/xl ;

    s->uMatrix.daa[2][2] = t;
    s->uMatrix.daa[k-1][k-1] = t;
    s->uMatrix.daa[2][k-1] = -t;
    s->uMatrix.daa[k-1][2] = -t;

    t = gj / xl;

    s->uMatrix.daa[3][3] = s->uMatrix.daa[ii-1][ii-1] =  t;
    s->uMatrix.daa[3][ii-1] = s->uMatrix.daa[ii-1][3] = -t;

    t = (eiy+ eiy) / xl;

    s->uMatrix.daa[4][4] = s->uMatrix.daa[jj-1][jj-1] = t+t;
    s->uMatrix.daa[4][jj-1] = s->uMatrix.daa[jj-1][4] = t;

    t = 6 * eiy/xl/xl;

    s->uMatrix.daa[2][4] = s->uMatrix.daa[4][2] = -t;
    s->uMatrix.daa[2][jj-1] = s->uMatrix.daa[jj-1][2] = -t;
    s->uMatrix.daa[k-1][4] = s->uMatrix.daa[4][k-1] =  t;
    s->uMatrix.daa[jj-1][k-1] = s->uMatrix.daa[k-1][jj-1] =  t;

    t = (eiz+ eiz) / xl;

    s->uMatrix.daa[5][5] = s->uMatrix.daa[kk-1][kk-1] = t+t;
    s->uMatrix.daa[5][kk-1] = s->uMatrix.daa[kk-1][5] = t;

    t = 6 * eiz/xl/xl;

    s->uMatrix.daa[1][5] = s->uMatrix.daa[5][1] = t;
    s->uMatrix.daa[1][kk-1] = s->uMatrix.daa[kk-1][1] = t;
    s->uMatrix.daa[5][j-1] = s->uMatrix.daa[j-1][5] = -t;
    s->uMatrix.daa[j-1][kk-1] = s->uMatrix.daa[kk-1][j-1] = -t;

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

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

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

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

       UnitsCopy( &(s->spColUnits[6]), &(s->spColUnits[0]) ); 
       UnitsCopy( &(s->spColUnits[7]), &(s->spColUnits[1]) ); 
       UnitsCopy( &(s->spColUnits[8]), &(s->spColUnits[2]) ); 
       UnitsCopy( &(s->spColUnits[9]), &(s->spColUnits[3]) ); 
       UnitsCopy( &(s->spColUnits[10]), &(s->spColUnits[4]) ); 
       UnitsCopy( &(s->spColUnits[11]), &(s->spColUnits[5]) ); 

       UnitsCopy( &(s->spRowUnits[6]), &(s->spRowUnits[0]) ); 
       UnitsCopy( &(s->spRowUnits[7]), &(s->spRowUnits[1]) ); 
       UnitsCopy( &(s->spRowUnits[8]), &(s->spRowUnits[2]) ); 
       UnitsCopy( &(s->spRowUnits[9]), &(s->spRowUnits[3]) ); 
       UnitsCopy( &(s->spRowUnits[10]), &(s->spRowUnits[4]) ); 
       UnitsCopy( &(s->spRowUnits[11]), &(s->spRowUnits[5]) ); 

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

    /* ===================== */
    /* Rotate stiff matrix   */
    /* ===================== */

    s->uMatrix.daa = (double **) rotate3d(s->uMatrix.daa, rot, no_dof);

#ifdef DEBUG
  printf("flag: in beamst3d() : STIFF after rotation \n");
  MatrixPrintIndirectDouble( s );
#endif
    
     return(s);
}


/* =============================================================== */
/* function   beamms3d()                                           */
/* 3D Frame Mass Matrix: This function returns a 12x12 rotated     */
/* mass matrix for 3D frame elmts.                                 */
/*                                                                 */
/* Input:                                                          */
/*   mbar->Density * X_area ;                                      */ 
/*   s->matrix pointer;                                            */
/*   xl->length       ; rot->rotation_matrix(12x12)                */
/*   nst->num_of_deg_of free per element; ndf->dof's per node;     */
/* =============================================================== */

#ifdef __STDC__
MATRIX *beamms3d(ARRAY *p, MATRIX *s, int mtype, double mbar, double xl, double rg , double **rot, int nst, int no_dof)
#else
MATRIX *beamms3d(p, s, mtype, mbar, xl, rg , rot, nst, no_dof)
ARRAY  *p;
MATRIX *s;
double **rot;
double mbar, rg , xl;
int    nst, no_dof, mtype;
#endif
{
int    n,m, i, j, k;
double t;
DIMENSIONS *d1,*d2,*d3;

#ifdef DEBUG
       printf("INFO >> In beamms3d() : no_dof= %d   \n",no_dof);
#endif

    switch(mtype) {
        case LUMPED:
	     t = mbar * xl/2;
             s->uMatrix.daa[0][0] = s->uMatrix.daa[1][1] = s->uMatrix.daa[2][2]  =  t;
	     s->uMatrix.daa[3][3] = t*rg*rg;
	     s->uMatrix.daa[4][4] = s->uMatrix.daa[5][5] = t*xl*xl/12.0;

             s->uMatrix.daa[no_dof][no_dof] = s->uMatrix.daa[no_dof+1][no_dof+1] = s->uMatrix.daa[no_dof+2][no_dof+2]  =  t;
	     s->uMatrix.daa[no_dof+3][no_dof+3] = t*rg*rg;
	     s->uMatrix.daa[no_dof+4][no_dof+4] = s->uMatrix.daa[no_dof+5][no_dof+5] = t*xl*xl/12.0;

             break;

        case CONSISTENT:

	    /* -------------------------------------- */
	    /* M terms =|  Mjj   Mjk |                */
	    /*          |  Mkj   Mkk |                */
	    /* -------------------------------------- */
	
	    /* -------------------------------------- */
	    /* Mjj  terms                             */
	    /* -------------------------------------- */

	    t = mbar*xl/420; 

	    s->uMatrix.daa[0][0]           =  t* 140;
	    s->uMatrix.daa[1][1] = s->uMatrix.daa[2][2] =  t* 156;
	    s->uMatrix.daa[3][3]           =  t* 140* rg * rg ;
	    s->uMatrix.daa[4][4] = s->uMatrix.daa[5][5] =  t* 4 * xl * xl;

	    s->uMatrix.daa[4][2] = s->uMatrix.daa[2][4] = -  t* 22 * xl;
	    s->uMatrix.daa[5][1] = s->uMatrix.daa[1][5] =    t* 22 * xl;

	    for(m=1; m <= no_dof; m++){
	        for(n = 1; n <= no_dof; n++){