fe_matrix.c

有限元分析源代码

   /* [c] : Assemble internal load vector from element level loads */

   for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) { 

   ep = &frame->element[elmt_no-1];
   if(displ_incr == (MATRIX *) NULL)
      ep->esp->state = 0; /* Calculate internal load as elastic behaviour */
   else
      ep->esp->state = 1; /* Need check for ELASTIC PLASTIC STATE */

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

   /*  Transfer fixed displacements */

   for(i=1; i <= array->nodes_per_elmt; i++) {

       k = 1;
       node_no = ep->node_connect[i-1];
       for(j = 1; j <= frame->no_dof; j++) {

           switch((int) array->nodes_per_elmt) {
               case 2: case 3:
                    jj = frame->node[node_no - 1].bound_id[j-1];  
                    if(jj > 0) {

                       /* displacement  */

                       array->displ->uMatrix.daa[j-1][i-1] = displ->uMatrix.daa[jj-1][0];
                       if( UNITS_SWITCH == ON ) {
                           UnitsCopy(&(array->displ->spRowUnits[j-1]), &(displ->spRowUnits[jj-1]));
                           UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
                       }

                       /* displacement increments */

                       if(displ_incr != (MATRIX *) NULL) {
                          array->displ_incr->uMatrix.daa[j-1][i-1] = displ_incr->uMatrix.daa[jj-1][0];
                          if(UNITS_SWITCH == ON) {
                             UnitsCopy(&(array->displ_incr->spRowUnits[j-1]), &(displ->spRowUnits[jj-1]));
                             UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
                          }
                       }
                    } else {

                       array->displ->uMatrix.daa[j-1][i-1] = frame->node[node_no -1].disp[j-1].value;

                       if(UNITS_SWITCH == ON ) {

                          /* displacement      */

                          UnitsCopy(&(array->displ->spRowUnits[j-1]), frame->node[node_no -1].disp[j-1].dimen);
                          UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));

                          /* displcement increments */

                          UnitsCopy(&(array->displ_incr->spRowUnits[j-1]), frame->node[node_no -1].disp[j-1].dimen);
                          UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
                       }
                    }
                    break;
               case 4: case 8:
                    jj = frame->node[node_no - 1].bound_id[j-1];
                    if(jj > 0) {

                       /* displacement  */

                       array->displ->uMatrix.daa[k-1][i-1] = displ->uMatrix.daa[jj-1][0];

                       if( UNITS_SWITCH == ON ) {
                           UnitsCopy(&(array->displ->spRowUnits[k-1]), &(displ->spRowUnits[jj-1]));
                           UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));
                       }

                       /* incremental displacement  */

                       if(displ_incr != (MATRIX *) NULL) {

                          array->displ_incr->uMatrix.daa[k-1][i-1] = displ_incr->uMatrix.daa[jj-1][0];

                          if(UNITS_SWITCH == ON) {
                             UnitsCopy(&(array->displ_incr->spRowUnits[k-1]), &(displ->spRowUnits[jj-1]));
                             UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
                          }
                       }
                    } else {
                       array->displ->uMatrix.daa[k-1][i-1] = frame->node[node_no -1].disp[j-1].value;

                       if(UNITS_SWITCH == ON ) {

                       /* displacement    */

                       UnitsCopy(&(array->displ->spRowUnits[k-1]), frame->node[node_no -1].disp[j-1].dimen);
                       UnitsCopy(&(array->displ->spColUnits[i-1]), &(displ->spColUnits[0]));

                       /* displcement increments */

                       UnitsCopy(&(array->displ_incr->spRowUnits[k-1]), frame->node[node_no-1].disp[j-1].dimen);
                       UnitsCopy(&(array->displ_incr->spColUnits[i-1]), &(displ->spColUnits[0]));
                       }
                    }
                    k = k + 1;
                    break;
               default:
                    break;
            }
       }
       }

       array = elmlib(array, LOAD_MATRIX);
       Ld    = Destination_Array(frame, elmt_no);



       if( UNITS_SWITCH == ON ) 
           ZeroUnits(&(load->spColUnits[0]));

       for(j = 1; j <= Ld[0]; j++) {
           k = Ld[j];
           if(k > 0) {
              load->uMatrix.daa[k-1][0] += array->nodal_loads[j-1].value;
              if( UNITS_SWITCH == ON )
                  UnitsCopy(&(load->spRowUnits[k-1]), array->nodal_loads[j-1].dimen);
           }
       }

#ifdef DEBUG
       printf("DIAGNOSTIC : elmt node %4d\n", elmt_no );
       dMatrixPrint("elmt diplacements", array->displ->uMatrix.daa, 3,2);
       printf("destination      :");
       for(itemp = 1; itemp <= Ld[0]; itemp++) {
           printf(" %10d : ", Ld[itemp]);
       }
       printf("\n");
       printf("elmt nodal loads :");
       for(itemp = 1; itemp <= Ld[0]; itemp++) {
           printf(" %10.2e : ", array->nodal_loads[itemp-1].value);
       }
       printf("\n");
       dMatrixPrint("Internal Load", load->uMatrix.daa, TNEQ, 1);
#endif

       free((char *) Ld);
   }

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

    return(load);
}


/* 
 *  =================================================
 *  UTILITY Functions for Stiffness and Mass matrices
 *  =================================================
 */ 

/* 
 *  -----------------------------------------
 *  [a] : Transfer Destination Array to *Ld[]
 *  -----------------------------------------
 */ 

#ifdef __STDC__
int *Destination_Array(EFRAME *frame, int elmt_no)
#else
int *Destination_Array(frame, elmt_no)
EFRAME *frame;
int    elmt_no;
#endif
{
ELEMENT *el;
int     i,no_dof_per_elmt, *Ld;
      
#ifdef DEBUG
    printf(" enter Destination_Array() \n");
    printf(" in Destination_Array() : \n");
    printf("                        : frame->no_dof = %d \n", frame->no_dof);
    printf("                        : frame->no_nodes_per_elmt = %d \n", frame->no_nodes_per_elmt);
#endif

    no_dof_per_elmt  = frame->no_dof*frame->no_nodes_per_elmt;

#ifdef DEBUG
    printf("                        : no_dof_per_elmt = %d \n", no_dof_per_elmt);
#endif

    Ld = iVectorAlloc((no_dof_per_elmt + 1));

    el = &frame->element[elmt_no - 1];
    for(i = 1; i <= no_dof_per_elmt + 1; i++) {
         Ld[i-1] = el->d_array[i-1];
    }

#ifdef DEBUG
    printf(" leaving  Destination_Array()  \n");
#endif
    return(Ld);
}

/* 
 *  ----------------------------------------------------
 *  [c] : Compute Destination Array *Ld[] for Rigid Body
 *  ----------------------------------------------------
 */ 

#ifdef __STDC__
int *Destination_Array_for_Rigid_Body(EFRAME *frp, ARRAY *p, int elmt_no, int dflg)
#else
int *Destination_Array_for_Rigid_Body(frp, p, elmt_no, dflg)
EFRAME       *frp;
ARRAY          *p;
int elmt_no, dflg;
#endif
{
ELEMENT                        *el;
int      no_nodes_per_elmt, no_dof;
int i,j,k,jj,ii,iid,mate_no, rb_ty;
int                           *vld;

       no_nodes_per_elmt = p->nodes_per_elmt;                  /* MIN(p->nodes_per_elmt, frp->no_nodes_per_elmt);*/
       no_dof            = p->dof_per_node;                    /* MIN(p->dof_per_node,frp->no_dof);*/
       vld     = iVectorAlloc(no_dof*no_nodes_per_elmt);

       rb_ty = frp->rigid[elmt_no -1].rb_type ;
       ii    = frp->rigid[elmt_no -1].in[1] ; /* node no of node conn by rbody for dofs */

       for(j=1;j<=no_dof;j++){     /* loop over dofs */
           jj = frp->rigid[elmt_no -1].rest_dof[j];
           if(jj >= 1 && rb_ty != 7) /* dofs not connected */
              vld[j] = 0;           /*k;*/
           else { /* jj=0 i:e; rest dofs & for  rb = 7; storey type rbody */
              k  = frp->node[ii - 1].bound_id[j]; /* equation nos assigned to dofs */
              vld[j] = k;
           }
       } 

       vld[0] = no_dof*no_nodes_per_elmt;

    return(vld);
}


/* 
 *  ------------------------------------------------------------------
 *  Assemble Global  Matrix : Stiffness and Mass  
 *  gmatrix : Global matrix  TNEQxTNEQ
 *  Ge      : Local matrix   array->size_of_stiffxarray->size_of_stiff
 *  ------------------------------------------------------------------
 */ 

#ifdef __STDC__
MATRIX *Assemble_Global(MATRIX *gmatrix, EFRAME *frame, int *Ld, MATRIX *Ge)
#else
MATRIX *Assemble_Global(gmatrix, frame, Ld, Ge)
MATRIX *gmatrix;
EFRAME *frame;
int    *Ld;
MATRIX *Ge;
#endif
{
ELEMENT          *el;
int j,l,k,m, size_ld;
int length1, length2;
int       iMin, iMax;
int     UNITS_SWITCH;
int linked, count, ik;

    size_ld = Ld[0];
    UNITS_SWITCH = CheckUnits();

#ifdef DEBUG
       printf("*** Enter Assemble_Global() : length of Ld = %3d : ", size_ld);
       for(j = 1; j <= size_ld; j++) {
           printf(" %3d ", Ld[j]);
       }
       printf("\n");
#endif

    for(j = 1; j <= size_ld; j++) {
        k = Ld[j];

        count = 0;
        for(ik = 1; ik <= size_ld; ik++) {
            if(k == Ld[ik]) count = count + 1;
        }

        linked = (count > 1);
        if(k > 0 && linked == 0) {
           for(l = 1; l <= size_ld; l++){
               m = Ld[l];
               if(m > 0){
	          switch( gmatrix->eRep ) {
	            case  INDIRECT :
                        gmatrix->uMatrix.daa[k-1][m-1] += Ge->uMatrix.daa[j-1][l-1];
	                break;
	            case  SKYLINE :
	                if( l >= j ) {
	                    iMin = (int) MIN( k, m );
	                    iMax = (int) MAX( k, m );
	                    if( (iMax-iMin+1) <= gmatrix->uMatrix.daa[iMax-1][0] )
                                gmatrix->uMatrix.daa[iMax-1][iMax-iMin+1]
                                += Ge->uMatrix.daa[j-1][l-1];
	                    else
	                        gmatrix->uMatrix.daa[iMax-1][iMax-iMin+1]
                                =  Ge->uMatrix.daa[j-1][l-1];
	                }
	                break;
	            default:
                        FatalError("In Assemble_Global() : Undefined spA->eRep",
                                    (char *) NULL);
	                break;
	          }
                  
                  if( UNITS_SWITCH == ON ) {
                     UnitsCopy( &(gmatrix->spColUnits[m-1]), &(Ge->spColUnits[l-1]) );
                     UnitsCopy( &(gmatrix->spRowUnits[k-1]), &(Ge->spRowUnits[j-1]) );
                  }
               }
           }
#ifdef DEBUG
           printf(" local dof : j = %d \t global dof : k = %d \n", j,k); 
#endif
        }
    }

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

    return(gmatrix);
}



/* 
 *  -----------------------------------------------
 *  Assemble Global Load: Equiv Nodal Load  
 *  -----------------------------------------------
 */ 

#ifdef __STDC__
MATRIX *Assemble_Global_Load(MATRIX *gmatrix, EFRAME *frame, int *Ld, MATRIX *Ge)
#else
MATRIX *Assemble_Global_Load(gmatrix, frame, Ld, Ge)
MATRIX *gmatrix;
EFRAME *frame;
int    *Ld;
MATRIX *Ge;
#endif
{
ELEMENT          *el;
int j,l,k,m, size_ld;
int           length;
int     UNITS_SWITCH;

    size_ld = Ld[0];
    UNITS_SWITCH = CheckUnits();

#ifdef DEBUG
       printf("*** Enter Assemble_Global_Load() : length of Ld = %3d : ", size_ld);
#endif

    for(j = 1; j <= size_ld; j++) {
        k = Ld[j];
        if(k > 0) {
          gmatrix->uMatrix.daa[k-1][0] += Ge->uMatrix.daa[j-1][0];
          if( UNITS_SWITCH == ON ) {
              UnitsCopy( &(gmatrix->spRowUnits[k-1]),  &(Ge->spRowUnits[j-1]) );
          }
        }
    }
    if( UNITS_SWITCH == ON )
       ZeroUnits( &(gmatrix->spColUnits[0]) );

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

    return(gmatrix);
}

/* ======================================================= */
/* Set Element Matrix - Stiffness/Mass for Feap problem    */
/* Input  - frame pointer, element number,p array          */
/* Output - Local Element Stiffness/Mass                   */
/* ======================================================= */

#ifdef __STDC__
MATRIX *Element_Matrix(ARRAY *p, int isw)
#else
MATRIX *Element_Matrix(p, isw)
ARRAY    *p;
int     isw;
#endif
{
   p = elmlib(p,isw);
   return(p->stiff);
}

/* ========================================================== */
/* Set Element Equivalent Load -                              */
/* Input  - p working array , task id :isw = EQUIV_NODAL_LOAD */
/* Output - Local Element Equivalent nodal load               */
/* ========================================================== */

#ifdef __STDC__
MATRIX *Element_Equiv(ARRAY *p, int isw)
#else
MATRIX *Element_Equiv(p, isw)
ARRAY    *p;
int     isw;
#endif
{