elmt_fiber3d.c

有限元程序

             /* put displacement matrix to one column form */
             /* p->displ=[px1,px2; py1,py2; pz1,pz2; rx1,rx2; ry1,ry2; rz1,rz2]6x2    */
             /* => temp_m1=pt=[px1;py1;pz1;rx1;ry1;rz1; px2;py2;pz2;rx2;ry2;rz2]12x1 */
             temp_m1  = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 12, 1 );
             for( ii=0 ; ii < p->dof_per_node ; ++ii )
                for( jj=0 ; jj < p->nodes_per_elmt ; ++jj )
                   temp_m1->uMatrix.daa[ii+jj*p->dof_per_node][0] = p->displ->uMatrix.daa[ii][jj];

             L = Element_Transformation_3d( p->coord, elmt_length );
             q = MatrixMult( L, temp_m1 );
             Q = MatrixMult( K, q );
             MatrixFree( temp_m1 );
          }
          else if( p->elmt_state == 1 ) {
             L = Element_Transformation_3d( p->coord, elmt_length );
             Q = p->Q_saved;
          }

          /* Calculate torsional stiffness */
          /* Assuming plane remain plane, so twisting stiffness is independent to others */
          if ( (G=p->work_material[1].value) <= 0.0 )
             G = p->work_material[0].value/2.0/(1+p->work_material[4].value);
          J = p->work_section[12].value;
          bf = p->work_section[7].value;
          depth = p->work_section[9].value;

          /* 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 = (1.0-0.63*bf/depth)*bf*bf*bf*depth/3.0;
             else
                J = (1.0-0.63*depth/bf)*depth*depth*depth*bf/3.0;
          }

          /* Transform local coordinate system to global coordinate system */

          cx = (p->coord[0][1].value - p->coord[0][0].value)/elmt_length.value;
          cy = (p->coord[1][1].value - p->coord[1][0].value)/elmt_length.value;
          cz = (p->coord[2][1].value - p->coord[2][0].value)/elmt_length.value;

          /* independent twisting stiffness */
          /* and element twising force in local coordinate */

          temp1 = G*J/elmt_length.value;
          temp2 = temp1*( cx*(p->displ->uMatrix.daa[3][1]-p->displ->uMatrix.daa[3][0]) + cy*(p->displ->uMatrix.daa[4][1]-p->displ->uMatrix.daa[4][0]) + cz*(p->displ->uMatrix.daa[5][1]-p->displ->uMatrix.daa[5][0]) );

          if( isw == LOAD_MATRIX ) {
             Ltrans = MatrixTranspose( L );
             temp_m1 = MatrixMult( Ltrans, Q );   /* element nodal force in global coord. */

             /* assign twisting force */
             temp_m1->uMatrix.daa[3][0]  = temp_m1->uMatrix.daa[3][0]  - cx*temp2;
             temp_m1->uMatrix.daa[4][0]  = temp_m1->uMatrix.daa[4][0]  - cy*temp2;
             temp_m1->uMatrix.daa[5][0]  = temp_m1->uMatrix.daa[5][0]  - cz*temp2;
             temp_m1->uMatrix.daa[9][0]  = temp_m1->uMatrix.daa[9][0]  + cx*temp2;
             temp_m1->uMatrix.daa[10][0] = temp_m1->uMatrix.daa[10][0] + cy*temp2;
             temp_m1->uMatrix.daa[11][0] = temp_m1->uMatrix.daa[11][0] + cz*temp2;
             MatrixFree( Ltrans );
          }
          if( isw == STRESS ) {
             R = Rigid_Body_Rotation_3d( elmt_length );
             Rtrans = MatrixTranspose( R );
             temp_m1 = MatrixMult( Rtrans, Q );   /* element nodal force in local coord. */

             /* assign twisting force */
             temp_m1->uMatrix.daa[3][0] = -temp2;
             temp_m1->uMatrix.daa[9][0] =  temp2;
             MatrixFree( R );
             MatrixFree( Rtrans );
          }
          MatrixFree( L );

          /* Assign force values */
          for( ii=0 ; ii < p->size_of_stiff ; ++ii )
             p->nodal_loads[ii].value = temp_m1->uMatrix.daa[ii][0];

          if( p->elmt_state == 0 ) {
             MatrixFree( K );
             MatrixFree( Q );
             MatrixFree( q );
          }
          MatrixFree( temp_m1 );

          /* Assign force units */
          if( UNITS_SWITCH == ON ) {
             SetUnitsOn();
             switch( UnitsType ) {
                case SI:
                case SI_US:
                   dp_force  = DefaultUnits("N");
                   dp_length = DefaultUnits("m");
                   break;
                case US:
                   dp_force  = DefaultUnits("lbf");
                   dp_length = DefaultUnits("in");
                   break;
             }
             dimen = UnitsMult( dp_force, dp_length );

             for( ii=1 ; ii <= 3 ; ++ii ) {
                UnitsCopy( p->nodal_loads[ii-1].dimen, dp_force );
                UnitsCopy( p->nodal_loads[ii-1+p->dof_per_node].dimen, dp_force );
                UnitsCopy( p->nodal_loads[ii-1+3].dimen, dimen );
                UnitsCopy( p->nodal_loads[ii-1+3+p->dof_per_node].dimen, dimen );
             }
          }

          if(isw == STRESS ) {
             xx  = 0.5*(p->coord[0][0].value + p->coord[0][1].value);   /* xx = 0.5(x1+x2) */
             yy  = 0.5*(p->coord[1][0].value + p->coord[1][1].value);   /* yy = 0.5(y1+y2) */
             zz  = 0.5*(p->coord[2][0].value + p->coord[2][1].value);   /* zz = 0.5(z1+z2) */
             printf("\n");
             printf("Elmt No %3d : ", p->elmt_no);
             switch( UNITS_SWITCH ) {
                case ON:
                   printf("Coords (X,Y,Z)= (%8.3f %s,%8.3f %s,%8.3f %s)\n\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);
                   /* 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/dimen->scale_factor, dimen->units_name,
                         p->nodal_loads[4].value/dimen->scale_factor, dimen->units_name,
                         p->nodal_loads[5].value/dimen->scale_factor, dimen->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/dimen->scale_factor, dimen->units_name,
                         p->nodal_loads[10].value/dimen->scale_factor, dimen->units_name,
                         p->nodal_loads[11].value/dimen->scale_factor, dimen->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"); 

                   free((char *) dp_force->units_name);
                   free((char *) dp_length->units_name);
                   free((char *) dimen->units_name);
                   free((char *) dp_force);
                   free((char *) dp_length);
                   free((char *) dimen);
                   break;
                case OFF:
                   printf("Coords (X,Y,Z)= (%8.3f,%8.3f,%8.3f)\n\n", xx, yy, zz);
                   /* 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;
             }
            }
            break;  /* end of case STRESS, LOAD_MATRIX */

       case STRESS_MATRIX:

            /* save element nodal forces in working array */

            printf("--------------------------------------------------- \n");
            printf("In elmt_fiber3d.c : case STRESS_MATRIX              \n");
            printf("--------------------------------------------------- \n");
            printf("*** This block of code has not yet been implemented \n");
            printf("*** See elmt_fiber2d.c for details in 2-d case      \n");
            printf("*** Terminating Program Execution                   \n");
            exit(1);

            break;

       case MASS_MATRIX:
          if( p->work_section[6].value != 0.0 )    /* mbar = weight/gravity */
             mbar = p->work_section[6].value/9.80665;
          /* mbar = density * area */
          else if( (p->work_material[5].value > 0.0) && (p->work_section[10].value > 0.0) )
             mbar = p->work_material[5].value * p->work_section[10].value;
          else
             FatalError("\nError in input: Need density value to calculate mass matrix\n",(char *)NULL);

          cx = p->coord[0][1].value - p->coord[0][0].value;           /* Cos Term */
          cy = p->coord[1][1].value - p->coord[1][0].value;           /* Sin Term */
          cz = p->coord[2][1].value - p->coord[2][0].value;           /* Tan Term */
          p->length.value = elmt_length.value;

          /* 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 radius of gyration , rT  --  m  ,  in   */
          /* original version      :  rT = p->length.value/ 1.414; */
          bf = p->work_section[7].value;
          depth = p->work_section[9].value;
          J = p->work_section[12].value;
          rT = p->work_section[13].value;

          /* 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 = (1.0-0.63*bf/depth)*bf*bf*bf*depth/3.0;
             else
                J = (1.0-0.63*depth/bf)*depth*depth*depth*bf/3.0;
          }
          if( rT==0 && p->work_section[10].value!=0 )
             rT = sqrt( J / p->work_section[10].value );

          switch( p->type ) {  /* mass type */
             case LUMPED :
                /* use lumped mass matrix of FRAME_3D element */
                p->stiff = beamms3d(p, p->stiff,p->type, mbar, elmt_length.value, rT, rot, p->size_of_stiff, p->dof_per_node);
                MatrixFreeIndirectDouble(rot, p->size_of_stiff);
                break;

             case CONSISTENT :
             default :
                FatalError("FIBER_3D element only support LUMPED mass", (char *)NULL);
                break;
          }

          break;

       default:
          break;
    }

    return(p);
}

/*
 *  ==================================
 *  Print Fiber Element Properties
 *  ==================================
 */

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

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

     printf("             ");
     printf(" General Section and Material Properties\n");
     switch(UNITS_SWITCH) {
       case ON:
        UnitsSimplify( eap->work_material[0].dimen );
        UnitsSimplify( eap->work_material[1].dimen );
        UnitsSimplify( eap->work_material[2].dimen );
        UnitsSimplify( eap->work_material[3].dimen );
        UnitsSimplify( eap->work_material[5].dimen );
        UnitsSimplify( eap->work_material[10].dimen );
        UnitsSimplify( eap->work_material[11].dimen );
        UnitsSimplify( eap->work_section[10].dimen );
	if( eap->work_material[5].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Density                 =      %11.3e %s\n",
                           eap->work_material[5].value/eap->work_material[5].dimen->scale_factor,
                           eap->work_material[5].dimen->units_name);
	}