fe_matrix.c

有限元分析源代码

       printf("*** Leave Get_Node()\n");
#endif

    return(connect);
}

#ifdef __STDC__
MATRIX *Get_Displ(MATRIX *m1, MATRIX *m2)
#else
MATRIX *Get_Displ(m1,m2)
MATRIX *m1, *m2;
#endif
{
MATRIX *displ;
int    nodeno;
int     ii,jj;
int UNITS_SWITCH;

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

       nodeno = (int) m1->uMatrix.daa[0][0];
       displ  = MatrixAllocIndirect("Node Displ", DOUBLE_ARRAY, 1, frame->no_dof);
       UNITS_SWITCH = CheckUnits();

          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]), UNITS_TYPE);
                         }
                         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
                         RadUnitsSimplify( frame->node[nodeno-1].disp[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);
}

#ifdef __STDC__
MATRIX *Get_Stress(MATRIX *m1, MATRIX *m2)
#else
MATRIX *Get_Stress(m1,m2)
MATRIX *m1, *m2;
#endif
{
MATRIX *stress;
int    elmt_no;
int UNITS_SWITCH;

ELEMENT        *ep;
ELEMENT_ATTR   *eap;
int node_no, elmt_attr_no;
int  i,j,k,ii,jj;

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

       elmt_no = (int) m1->uMatrix.daa[0][0];
       stress  = MatrixAllocIndirect("Element Stress", DOUBLE_ARRAY, frame->no_nodes_per_elmt, frame->no_dof);
       UNITS_SWITCH = CheckUnits();

#ifdef DEBUG
        printf(" elmt No = %d \n", elmt_no);
#endif

           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;
                   }
               }
           }
	   PRINT_STRESS = OFF;
           array = elmlib(array, STRESS);

           for( i=1 ; i <= frame->no_dof ; i++ )
              for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
                  frame->element[elmt_no-1].rp->Forces->uMatrix.daa[i-1][j-1]
                  = array->nodal_loads[frame->no_dof*(j-1)+i-1].value;
           if( UNITS_SWITCH == ON ) {
              for( i=1 ; i <= frame->no_dof ; i++ )
                  UnitsCopy( &(frame->element[elmt_no-1].rp->Forces->spRowUnits[i-1]) , array->nodal_loads[i-1].dimen );
              for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
                  ZeroUnits( &(frame->element[elmt_no-1].rp->Forces->spColUnits[j-1]) );
           }

           for( i=1 ; i <= frame->no_dof ; i++ )
              for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
                  stress->uMatrix.daa[j-1][i-1] = frame->element[elmt_no-1].rp->Forces->uMatrix.daa[i-1][j-1];
           if( UNITS_SWITCH == ON ) {
              for( i=1 ; i <= frame->no_dof ; i++ )
                  UnitsCopy( &(stress->spColUnits[i-1]), &(frame->element[elmt_no-1].rp->Forces->spRowUnits[i-1]) );
              for( j=1 ; j <= frame->no_nodes_per_elmt ; j++ )
                  ZeroUnits( &(stress->spRowUnits[j-1]) );
           }

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

    return(stress);
}

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

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

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

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

    return(gdof);
}

#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
/*
        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
*/

       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->spRowUnits[14]) );
          if( sap->tw.dimen != (DIMENSIONS *)NULL )
             UnitsCopy( &(section->spRowUnits[15]), sap->tw.dimen );
          else
             ZeroUnits( &(section->spRowUnits[15]) );
       }

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

    return(section);
}

#ifdef __STDC__
MATRIX *Get_Material(MATRIX *m)
#else
MATRIX *Get_Material(m)
MATRIX *m;
#endif
{
MATRIX *material;
int       elmtno;
int elmt_attr_no;
int           ii;
char       *name;
MATERIAL_ATTR *map;

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

/*
        QUANTITY                   E;    [0] Young's modulus
        QUANTITY                   G;    [1] Shear modulus
        QUANTITY                  fy;    [2] Yield stress
        QUANTITY                  ET;    [3] Tangent Young's Modulus
        double                    nu;    [4] Poission's ratio
        QUANTITY             density;    [5] 
        QUANTITY                  fu;    [6] Ultimate stress
        QUANTITY      *alpha_thermal;    [7] thermal expansion coefficient
*/
       elmtno = (int) m->uMatrix.daa[0][0];
       elmt_attr_no = frame->element[elmtno-1].elmt_attr_no;
       name = frame->eattr[elmt_attr_no-1].material;
       map = lookup(name)->u.map;
       material  = MatrixAllocIndirect("Elmt Material", DOUBLE_ARRAY, 10, 1);

       material->uMatrix.daa[0][0] = map->E.value;
       material->uMatrix.daa[1][0] = map->G.value;
       material->uMatrix.daa[2][0] = map->fy.value;
       material->uMatrix.daa[3][0] = map->ET.value;
       material->uMatrix.daa[4][0] = map->nu;
       material->uMatrix.daa[5][0] = map->density.value;
       material->uMatrix.daa[6][0] = map->fu.value;
       material->uMatrix.daa[7][0] = map->alpha_thermal[0].value;
       material->uMatrix.daa[8][0] = map->alpha_thermal[1].value;
       material->uMatrix.daa[9][0] = map->alpha_thermal[2].value;
       if( CheckUnits() == ON ) {
          ZeroUnits( &(material->spColUnits[0]) );
          if( map->E.dimen != (DIMENSIONS *)NULL )
             UnitsCopy( &(material->spRowUnits[0]), map->E.dimen );
          else
             ZeroUnits( &(material->spRowUnits[0]) );
          if( map->G.dimen != (DIMENSIONS *)NULL )
             UnitsCopy( &(material->spRowUnits[1]), map->G.dimen );
          else
             ZeroUnits( &(material->spRowUni