fe_matrix.c

有限元程序

       displ  = MatrixAllocIndirect("Node Displ", DOUBLE_ARRAY, 1, frame->no_dof);
       UNITS_SWITCH = CheckUnits();
       UnitsType    = CheckUnitsType();

       switch( UNITS_SWITCH ) {
          case ON:
               for( ii=0 ; ii < frame->no_dof ; ii++ ) {
                    jj = frame->node[nodeno-1].bound_id[ii];
                    if(jj > 0) {
                       if(m2->spRowUnits[jj-1].units_name != NULL ) {
                          UnitsTypeConvert(&(m2->spRowUnits[jj-1]), UnitsType);
                       }
                       RadUnitsSimplify( &(m2->spRowUnits[jj-1]) );
                       frame->node[nodeno-1].disp[ii].value = m2->uMatrix.daa[jj-1][0];
                       UnitsCopy( frame->node[nodeno-1].disp[ii].dimen, &(m2->spRowUnits[jj-1]) );
                    }
                    else {
                       if ( ii < frame->no_dimen ) {
                            UnitsCopy( frame->node[nodeno-1].disp[ii].dimen, 
                                       frame->node[nodeno-1].coord[ii].dimen );
                       }
                    }
               }
               break;
          case OFF:
               for( ii=0 ; ii<frame->no_dof ; ii++){
                    jj = frame->node[nodeno-1].bound_id[ii];
                    if(jj > 0)
                       frame->node[nodeno-1].disp[ii].value = m2->uMatrix.daa[jj-1][0];
               }
               break;
          default:
               break;
       } 

       if( UNITS_SWITCH == ON ) {
          for( ii=0 ; ii<frame->no_dof ; ii++ ) {
             displ->uMatrix.daa[0][ii] = frame->node[nodeno-1].disp[ii].value;
             UnitsCopy( &(displ->spColUnits[ii]), frame->node[nodeno-1].disp[ii].dimen );
          }
          ZeroUnits( &(displ->spRowUnits[0]) );
       }
       else {
          for( ii=0 ; ii<frame->no_dof ; ii++ )
             displ->uMatrix.daa[0][ii] = frame->node[nodeno-1].disp[ii].value;
       }

#ifdef DEBUG
       printf("*** Leave Get_Displ()\n");
#endif

    return(displ);
}

/* 
 *  ================================================================= 
 *  Get_Stress() : Return stresses within an element.
 *
 *  Input  : Matrix *m      -- node no.
 *  Output : Matrix *stress -- matrix of stresses within element.
 *  ================================================================= 
 */ 

#ifdef __STDC__
MATRIX *Get_Stress(MATRIX *m1, MATRIX *m2)
#else
MATRIX *Get_Stress(m1,m2)
MATRIX *m1, *m2;
#endif
{
MATRIX        *stress;
ELEMENT           *ep;
ELEMENT_ATTR     *eap;
DIMENSIONS *dp_length;
int UNITS_SWITCH;
int      elmt_no;
int node_no, elmt_attr_no;
int  i,j,k,ii,jj, iNoCols;

    elmt_no = (int) m1->uMatrix.daa[0][0];

#ifdef DEBUG
       printf("*** Enter Get_Stress() : elmt_no = %d\n", elmt_no );
#endif

    /* Allocate working array for elmt_no and assign element properties */

    UNITS_SWITCH = CheckUnits();
    array = Assign_p_Array(frame, elmt_no, array, STRESS);
    array = elmlib(array, PROPTY);

    /* Transfer Fixed Displacements */

    ep           = &frame->element[elmt_no-1];
    elmt_attr_no = ep->elmt_attr_no;  
    eap          = &frame->eattr[elmt_attr_no-1];

    for(i=1; i <= array->nodes_per_elmt; i++) {
        k = 1; 
        node_no = ep->node_connect[i-1];
        for(j = 1; j <= array->dof_per_node; j++) {
            switch( (int) array->nodes_per_elmt ) {
                case 2:
                case 3:
                     ii = eap->map_ldof_to_gdof[j-1];
                     jj = frame->node[node_no - 1].bound_id[ii-1];
                     if(jj > 0) {
                        array->displ->uMatrix.daa[j-1][i-1] = m2->uMatrix.daa[jj-1][0];
                        if( UNITS_SWITCH == ON ) {
                            UnitsCopy(&(array->displ->spRowUnits[j-1]), &(m2->spRowUnits[jj-1]));
                            ZeroUnits(&(array->displ->spColUnits[i-1]));
                        }
                     } else {
                        array->displ->uMatrix.daa[j-1][i-1]
                           = frame->node[node_no -1].disp[ii-1].value;
                        if( UNITS_SWITCH == ON ) {
                            UnitsCopy(&(array->displ->spRowUnits[j-1]),
                                      frame->node[node_no -1].disp[ii-1].dimen);
                            ZeroUnits(&(array->displ->spColUnits[i-1]));
                        }
                     }
                     break;
                case 4:
                case 8:
                     ii = eap->map_ldof_to_gdof[k-1];
                     jj = frame->node[node_no - 1].bound_id[ii-1];
                     if(jj > 0) {
                        array->displ->uMatrix.daa[k-1][i-1] = m2->uMatrix.daa[jj-1][0];
                        if( UNITS_SWITCH == ON ) {
                            UnitsCopy(&(array->displ->spRowUnits[k-1]), &(m2->spRowUnits[jj-1]));
                            ZeroUnits( &(array->displ->spColUnits[i-1]) );
                        }
                     } else {
                        array->displ->uMatrix.daa[k-1][i-1]
                           = frame->node[node_no -1].disp[ii-1].value;
                        if( UNITS_SWITCH == ON ) {
                            UnitsCopy( &(array->displ->spRowUnits[k-1]),
                                       frame->node[node_no -1].disp[ii-1].dimen);
                            ZeroUnits( &(array->displ->spColUnits[i-1]) );
                        }
                     }
                     k = k + 1;
                     break;
                 default:
                     break;
           }
       }
   }

   /* Allocate memory for the array of element stresses/internal forces. */

   switch ( frame->no_dimen ) {
       case 1:
            iNoCols = 2;
            break;
       case 2:
            iNoCols = 5;
            break;
       case 3:
            iNoCols = 9;
            break;
       default:
            break;
   }

   stress = MatrixAllocIndirect("Element Stress", DOUBLE_ARRAY,
                                array->nodes_per_elmt, iNoCols );

   /* Retrieve element level stresses/forces */

#ifdef DEBUG
       printf("*** Go to elmlib( array, STRESS_MATRIX )\n");
#endif

   array = elmlib(array, STRESS_MATRIX );

#ifdef DEBUG
       printf("*** Back from elmlib( array, STRESS_MATRIX )\n");
#endif

   /* Transfer "units" to "stress" matrix */

   if( UNITS_SWITCH == ON ) {
       for( i = 1; i <= 4*frame->no_dimen - 3 ; i++ ) {
            UnitsCopy( &(stress->spColUnits[ i - 1 ]),
                       &(array->stress->spColUnits[i-1]) );
       }

       for( j = 1; j <= array->nodes_per_elmt ; j++ ) {
            ZeroUnits( &(stress->spRowUnits[j-1]) );
       }
   }

   /* Transfer "coordinate" and "force/stress values" to "stress" matrix */

   for( i = 1; i <= array->nodes_per_elmt; i++ )
   for( j = 1; j <= 4*frame->no_dimen - 3 ; j++ )
        stress->uMatrix.daa[i-1][j-1] = array->stress->uMatrix.daa[i-1][j-1];

#ifdef DEBUG
       printf("*** Leave Get_Stress()\n");
#endif

    return(stress);
}

/* 
 *  ================================================================= 
 *  Get_Dof() : Return dof associated with a particular node.
 *
 *  Input  : Matrix *m      -- node no.
 *  Output : Matrix *gdof   -- matrix of degrees of freedom.
 *  ================================================================= 
 */ 

#ifdef __STDC__
MATRIX *Get_Dof(MATRIX *m)
#else
MATRIX *Get_Dof(m)
MATRIX *m;
#endif
{
MATRIX  *gdof;
int    nodeno;
int        ii;

   nodeno = (int) m->uMatrix.daa[0][0];
   gdof   = MatrixAllocIndirect("Node Global Dof", DOUBLE_ARRAY, 1, frame->no_dof);

   if( CheckUnits() == ON ) {
       for( ii=0 ; ii<frame->no_dof ; ii++ ) {
            gdof->uMatrix.daa[0][ii] = (double) frame->node[nodeno-1].bound_id[ii];
            ZeroUnits( &(gdof->spColUnits[ii]) );
       }
       ZeroUnits( &(gdof->spRowUnits[0]) );
   } else {
       for( ii=0 ; ii<frame->no_dof ; ii++ )
             gdof->uMatrix.daa[0][ii] = (double) frame->node[nodeno-1].bound_id[ii];
   }

   return(gdof);
}


/* 
 *  ================================================================= 
 *  Get_Stiffness() : Return element stiffness matrix 
 *
 *  Input  : Matrix *m      -- element no.
 *  Output : Matrix *stiff  -- element stiffness matrix.
 *
 *  Note : This routine is a real hack (need to tidy up later).
 *  ================================================================= 
 */ 

#ifdef __STDC__
MATRIX *Get_Stiffness(MATRIX *m)
#else
MATRIX *Get_Stiffness(m)
MATRIX *m;
#endif
{
MATRIX      *stiff;
MATRIX         *Ke; 
MATRIX      *Ksave; 
int        elmt_no;
int        iElmtNo;

   iElmtNo = (int) m->uMatrix.daa[0][0];

   /* [a] : Compute element level stiffness matrices */

   array = Assign_p_Array(frame, iElmtNo, array, STIFF);
   array = Assign_p_Array(frame, iElmtNo, array, PROPTY);
   array = Element_Property(array); 

   /* [b] : Compute and save element level stiffness matrices */

   Ke    = Element_Matrix(array, STIFF);
   Ksave = MatrixCopy( Ke );

   return( Ksave );
}


/* 
 *  ==================================================================== 
 *  Get_Section() : Return section properties of a particular element.
 *
 *  Input  : Matrix *m       -- elementno.
 *  Output : Matrix *section -- matrix of section properties.
 *  ==================================================================== 
 *  The details are:
 *
 *      QUANTITY  Ixx;              [0]
 *      QUANTITY  Iyy;              [1]
 *      QUANTITY  Izz;              [2]
 *      QUANTITY  Ixz;              [3]
 *      QUANTITY  Ixy;              [4]
 *      QUANTITY  Iyz;              [5]
 *      QUANTITY  weight;           [6]  Section weight
 *      QUANTITY  bf;               [7]  Width of flange
 *      QUANTITY  tf;               [8]  thickness of flange
 *      QUANTITY  depth;            [9]  Section depth
 *      QUANTITY  area;             [10] Section area
 *      QUANTITY  plate_thickness;  [11]
 *      QUANTITY  tor_const;        [12] Torsional Constant J
 *      QUANTITY  rT;               [13] Section radius of gyration
 *      QUANTITY  width;            [14] Section width
 *      QUANTITY  tw;               [15] Thickness of web
 *      double    ks;               [16] Shear Section Correction Factor
 *  ==================================================================== 
 */ 

#ifdef __STDC__
MATRIX *Get_Section(MATRIX *m)
#else
MATRIX *Get_Section(m)
MATRIX *m;
#endif
{
MATRIX  *section;
int       elmtno;
int elmt_attr_no;
char       *name;
SECTION_ATTR *sap;

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

   elmtno       = (int) m->uMatrix.daa[0][0];
   elmt_attr_no = frame->element[elmtno-1].elmt_attr_no;
   name         = frame->eattr[elmt_attr_no-1].section;
   sap          = lookup(name)->u.sap;
   section      = MatrixAllocIndirect("Elmt Section", DOUBLE_ARRAY, 16, 1);

   section->uMatrix.daa[0][0]  = sap->Ixx.value;
   section->uMatrix.daa[1][0]  = sap->Iyy.value;
   section->uMatrix.daa[2][0]  = sap->Izz.value;
   section->uMatrix.daa[3][0]  = sap->Ixz.value;
   section->uMatrix.daa[4][0]  = sap->Ixy.value;
   section->uMatrix.daa[5][0]  = sap->Iyz.value;
   section->uMatrix.daa[6][0]  = sap->weight.value;
   section->uMatrix.daa[7][0]  = sap->bf.value;
   section->uMatrix.daa[8][0]  = sap->tf.value;
   section->uMatrix.daa[9][0]  = sap->depth.value;
   section->uMatrix.daa[10][0] = sap->area.value;
   section->uMatrix.daa[11][0] = sap->plate_thickness.value;
   section->uMatrix.daa[12][0] = sap->tor_const.value;
   section->uMatrix.daa[13][0] = sap->rT.value;
   section->uMatrix.daa[14][0] = sap->width.value;
   section->uMatrix.daa[15][0] = sap->tw.value;

   if( CheckUnits() == ON ) {
       ZeroUnits( &(section->spColUnits[0]) );
       if( sap->Ixx.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[0]), sap->Ixx.dimen );
       else
           ZeroUnits( &(section->spRowUnits[0]) );
       if( sap->Iyy.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[1]), sap->Iyy.dimen );
       else
           ZeroUnits( &(section->spRowUnits[1]) );
       if( sap->Izz.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[2]), sap->Izz.dimen );
       else
           ZeroUnits( &(section->spRowUnits[2]) );
       if( sap->Ixz.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[3]), sap->Ixz.dimen );
       else
           ZeroUnits( &(section->spRowUnits[3]) );
       if( sap->Ixy.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[4]), sap->Ixy.dimen );
       else
           ZeroUnits( &(section->spRowUnits[4]) );
       if( sap->Iyz.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[5]), sap->Iyz.dimen );
       else
           ZeroUnits( &(section->spRowUnits[5]) );
       if( sap->weight.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[6]), sap->weight.dimen );
       else
           ZeroUnits( &(section->spRowUnits[6]) );
       if( sap->bf.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[7]), sap->bf.dimen );
       else
           ZeroUnits( &(section->spRowUnits[7]) );
       if( sap->tf.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[8]), sap->tf.dimen );
       else
           ZeroUnits( &(section->spRowUnits[8]) );
       if( sap->depth.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[9]), sap->depth.dimen );
       else
           ZeroUnits( &(section->spRowUnits[9]) );
       if( sap->area.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[10]), sap->area.dimen );
       else
           ZeroUnits( &(section->spRowUnits[10]) );
       if( sap->plate_thickness.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[11]), sap->plate_thickness.dimen );
       else
           ZeroUnits( &(section->spRowUnits[11]) );
       if( sap->tor_const.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[12]), sap->tor_const.dimen );
       else
           ZeroUnits( &(section->spRowUnits[12]) );
       if( sap->rT.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[13]), sap->rT.dimen );
       else
           ZeroUnits( &(section->spRowUnits[13]) );
       if( sap->width.dimen != (DIMENSIONS *)NULL )
           UnitsCopy( &(section->spRowUnits[14]), sap->width.dimen );
       else
           ZeroUnits( &(section->spRowUnit