fe_matrix.c

有限元分析源代码

   p = elmlib(p,isw);
   return(p->equiv_nodal_load);
}

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

#ifdef __STDC__
QUANTITY *Element_Vector(ARRAY *p, int isw)
#else
QUANTITY *Element_Vector(p, isw)
ARRAY     *p;
int      isw;
#endif
{

  /*printf("*** In Element_VECTOR. Flag 1. Elem_# = %d\n",p->elmt_no); */
  p = elmlib(p,isw);
  return(p->nodal_loads);
}

/* ===================================================================== */
/* Rigidbody_conn                                                        */
/* Input - frame pointer frp, array ptr checks if any node of an element */
/*         is connected to a rigid body. It labels it by adding          */
/*         a value of (one thousand * rigid_body no) to such nodes       */
/*                                                                       */
/* NOTE : Nodes attached to rigid body are scaled by a factor 1,000,000  */
/*        in the working-element p->ix[] data structure.                 */
/* ===================================================================== */

#ifdef __STDC__
int Rigidbody_conn(EFRAME *frp, ARRAY *p)
#else
int Rigidbody_conn(frp,p)
EFRAME *frp;
ARRAY  *p;
#endif
{
int i,j,k,nn, no_node,out;
out = NO;

    for(i = 1; i<=frp->no_rigid;i++) {
        no_node = frp->rigid[i-1].nodes_conn;

        for(j = 1; j<=no_node;j++){
            nn = frp->rigid[i-1].in[j];
            for(k = 1; k<=p->nodes_per_elmt;k++) {
                if(p->node_connect[k] == nn){
                   p->node_connect[k] = p->node_connect[k] + i * 1000000;
                   frp->node[nn -1].rb_num = i;
                   out = YES;
                }
            }
        }
    }

    return(out);
}


/* =================================================================== */
/* Transform Stiff Matrix of member connected to Rigid Body            */
/* transforms action at j of member to working pointt f of Rigid Body. */
/*                                                                     */
/* Transforms element S. Matrix Ke to T*Ke*Tt                          */
/* Input - EFRAME *frp;  ARRAY *p;                                     */
/*         Member Stiff Matrix MS                                      */
/* Output- Modified Member Stiff Matrix MS                             */
/* =================================================================== */

#ifdef __STDC__
MATRIX *Transform_Stiff_Matrix(EFRAME *frp, ARRAY *p, MATRIX *MS)
#else
MATRIX *Transform_Stiff_Matrix(frp, p, MS)
EFRAME     *frp;
ARRAY        *p;
MATRIX      *MS;
#endif
{
int    nn,i,j,k,rb_no,node_no,rb_ty;
double **Tfj, **Ti, **TiT, **ktit;
double cgx,cgy,cgz;

    Tfj = (double **)MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);
    Ti  = (double **)MatrixAllocIndirectDouble(p->size_of_stiff,p->size_of_stiff);
    TiT = (double **)MatrixAllocIndirectDouble(p->size_of_stiff,p->size_of_stiff);

    for(i = 1; i<=p->nodes_per_elmt; i++) { 
        nn = p->node_connect[i];
        if(nn > 1000000) {
           node_no = nn - rb_no * 1000000;
           rb_no = nn/1000000; /* find rigid body number */

    /* ---- FIX LATER 
           if(frp->ndm == 3) {
               if(frp->basic_elmt_type == TRUSS_3D)
                 rb_ty = TRUSS_3D;
              else
                 rb_ty = FRAME_3D;
            }
           else 
                 rb_ty = frp->rigid[rb_no -1].rb_type; 
           */ 

           cgx = p->coord[1][i].value - frp->rigid[rb_no -1].xcg.value;
           cgy = p->coord[2][i].value - frp->rigid[rb_no -1].ycg.value;

/******** fix later
           if(frp->ndm >= 3)
              cgz = p->coord[3][i] - frp->rigid[rb_no -1].zcg;
	   else
              cgz = 0.0;
****************************/

           Tfj = (double **) Transformation_Matrix(Tfj,rb_ty, p->dof_per_node,cgx,cgy,cgz);
           Tfj = (double **)       Modify_T_Matrix(Tfj,  frp, p->dof_per_node, rb_no);
       }
       else 
         Tfj = Transformation_Matrix(Tfj,rb_ty, p->dof_per_node,0.0,0.0,0.0);

       /* get other terms by symm          */
       /* printf("node at end = %d\n", i); */

       for(j = 1; j<=p->dof_per_node; j++) {
          for(k = 1; k<=p->dof_per_node; k++) {
               Ti[(i-1)*p->dof_per_node+j][ (i-1)*p->dof_per_node+ k] = Tfj[j][k];
              TiT[(i-1)*p->dof_per_node+k][ (i-1)*p->dof_per_node+ j] = Tfj[j][k];
          }
       }
    }

   /* ------------------ */
   /* Transformed Matrix */
   /* ------------------ */

   ktit = (double **) dMatrixMult(MS->uMatrix.daa,p->size_of_stiff,p->size_of_stiff, TiT, p->size_of_stiff, p->size_of_stiff);
   MS->uMatrix.daa   = (double **) dMatrixMultRep(MS->uMatrix.daa, Ti, p->size_of_stiff, p->size_of_stiff, ktit, p->size_of_stiff, p->size_of_stiff);

   MatrixFreeIndirectDouble(Tfj, p->dof_per_node);
   MatrixFreeIndirectDouble(Ti, p->size_of_stiff);
   MatrixFreeIndirectDouble(TiT, p->size_of_stiff);
   MatrixFreeIndirectDouble(ktit, p->size_of_stiff);

   return(MS);
}


/* ============================================================ */
/* Transform Mass Matrix of RigidBody:                          */
/* Transforms R body Mass Matrix M =p->s to Mp = Tpc.M.Ttpc.    */
/* transforms action from                                       */
/*  mass centres   to  working point of Rigid Body              */
/*                                                              */
/*  here pcx,pcy,pcz are cordinates of mass centre c of RBODY   */
/*                       wrt working point p of RBODY           */
/* Input -       frame pointer frp, Array ptr p.                */
/*                Member Stiff Matrix MS                        */
/* Output-       Modified Member Stiff Matrix MS                */
/* ============================================================ */

#ifdef __STDC__
MATRIX *Transform_Rigid_Body_Mass_Matrix(EFRAME *frp, ARRAY *p, MATRIX *MS)
#else
MATRIX *Transform_Rigid_Body_Mass_Matrix(frp, p, MS)
EFRAME   *frp;
ARRAY    *p;
MATRIX   *MS;
#endif
{
int nn, i,j,k,rb_no,node_no,rb_ty;
double **Tpc, **Ttpc,**MTtpc;
double pcx,pcy,pcz;

     /*----------------------------------------------------------------------*/
     /* Alloc_Matrix                                                         */
     /*----------------------------------------------------------------------*/

        Tpc  = MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);
        Ttpc = MatrixAllocIndirectDouble(p->dof_per_node,p->dof_per_node);

     /*----------------------------------------------------------------------*/
     /* find rigid body number */
     /*----------------------------------------------------------------------*/

     rb_no = p->elmt_no;         
     rb_ty = frp->rigid[rb_no -1].rb_type;

     /*----------------------------------------------------------------------*/
     /* Coordinates of mass centre c  wrt working point p(geom centre)       */
     /*    = global coord of cg c - global coord of wkg pt  p                */ 
     /*----------------------------------------------------------------------*/

     pcx =     p->work_section[4].value - frp->rigid[rb_no -1].xcg.value;
     pcy =     p->work_section[5].value - frp->rigid[rb_no -1].ycg.value;
     pcz = 0.0;
/***********
     if( frp->ndm >= 3)
        pcz =    p->work_section[6] - frp->rigid[rb_no -1].zcg;
***************/
     Tpc = (double **) Transformation_Matrix(Tpc,rb_ty, p->dof_per_node,pcx,pcy,pcz);

     /*----------------------------------------------------------------------*/
     /* Modify depending on the dofs constrained   */
     /*----------------------------------------------------------------------*/

     Tpc = (double **) Modify_T_Matrix(Tpc, frp,p->dof_per_node,rb_no);

    /* Transpose of Tpc matrix */

     for(j = 1; j<=p->dof_per_node;j++) 
       for(k = 1; k<=p->dof_per_node;k++) 
          Ttpc[j][k] = Tpc[k][j];
      
      /* Product : M.Ttpc           */
      /* Product : Mp = Tpc. M.Ttpc */

      MTtpc = (double **) dMatrixMult(MS->uMatrix.daa, p->dof_per_node, p->dof_per_node, Ttpc, p->dof_per_node, p->dof_per_node);
      MS->uMatrix.daa    = (double **) dMatrixMultRep(MS->uMatrix.daa, Tpc, p->dof_per_node, p->dof_per_node, MTtpc, p->dof_per_node, p->dof_per_node);

      MatrixFreeIndirectDouble(Tpc, p->dof_per_node);
      MatrixFreeIndirectDouble(Ttpc, p->dof_per_node);
      MatrixFreeIndirectDouble(MTtpc, p->dof_per_node);

      return(MS);
}


/* ======================================================*/
/* Transformation matrix Tfj (=Tpj) for point j to f(=p) */ 
/* Action:  Ap = Tpj . Aj; Disp:    Dj = Ttpj . Dp       */
/*   -- for framed structures and conn rigid bodies      */
/* function double **Transformation_Matrix               */
/*                    (t,type,no_dof,cx,cy,cz)              */
/*  input;  t         :allocated matrix                  */
/*          type      : type of element                  */
/*          no_dof       : no of dof                        */
/*          cx,cy,cz  : coordinates of  point j wrt      */
/*                                    reference point p  */
/* =======================================================*/

#ifdef __STDC__
double **Transformation_Matrix(double **t, char *type, int no_dof, double cx, double cy, double cz)
#else
double **Transformation_Matrix(t, type,no_dof,cx,cy,cz)
double **t;
char   *type;
int    no_dof;
double cx,cy,cz;
#endif
{
int i,j,n;
     for(n = 0; n <= NO_ELEMENTS_IN_LIBRARY; n++) { 

         if(!strcmp(type, elmt_library[n].name))
         break;
       }

       for(i = 1; i <= elmt_library[n].no_dof; i++)
       for(j = 1; j <= elmt_library[n].no_dof; j++)
           t[i][j] = 0.0;

       if(elmt_library[n].no_dof == 3) {
        /*2d beam:2d bar:2d frame: grid*/

             t[1][1] = 1.0;
             t[2][2] = 1.0;
             if(!strcmp(type,"BEAM_2D"))
                t[2][1] = cx;
             else if(!strcmp(type,"TRUSS_2D")) {
                t[3][1] = -cy;
                t[3][2] =  cx;
             }
             else if(!strcmp(type,"FRAME_2D")){
                t[3][1] = -cy;
                t[3][2] =  cx;
                t[3][3] =  1;
             }
             else {                               /*  GRID_XY */
                t[1][3] =  cy;
                t[2][3] = -cx;
                t[3][3] =  1;
             }

          }

       if(elmt_library[n].no_dof == 6) {

             for(i = 1; i <= 6; i++)
                 t[i][i] = 1.0;
     
             if(!strcmp(type,"TABLE_XZ")) {   /* 3d space table o== r storey */
                t[3][1] =  cz;
                t[3][2] = -cx;
             }
             else if(!strcmp(type,"TRUSS_3D")) {   /* space truss */
                t[5][1] =  cz;
                t[6][1] = -cy;
                t[4][2] = -cz;
                t[6][2] =  cx;
                t[4][3] =  cy;
                t[5][3] = -cx;
             }
             else {                               /* space frames */
                for(i = 4; i <= 6; i++)
                    t[i][i] = 1.0;
                t[4][2] = -cz;
                t[4][3] =  cy;
                t[5][1] =  cz;
                t[5][3] = -cx;
                t[6][1] = -cy;
                t[6][2] =  cx;
             }
     }

     return(t);
}


/* 
 *  =====================================================================
 *  Modify Transformation_Matrix(T) : makes rows and columns to zero for
 *  dof not connected to rigid body. Diagonal terms are retained as unity 
 *  
 *  Input - frame pointer frp, prob case  
 *  =====================================================================
 */ 

#ifdef __STDC__
double **Modify_T_Matrix(double **T, EFRAME *frp, int no_dof, int rb_no)
#else
double **Modify_T_Matrix(T,frp, no_dof,rb_no)
double **T;
EFRAME *frp;
int    no_dof, rb_no;
#endif
{
int nn, i,j,k,ndof;
int rb_elmt_type;

    rb_elmt_type  = frp->rigid[rb_no -1].rb_type;

    /* For Storey type rigid body                     */
    /* all dofs(=3) assumed to be connected to R.Body */

    if(rb_elmt_type == 7) { /* Storey type rigid body ;dofs=3 */ 
                            /* all dofs are connected hence returned unaltered     */
                            /* printf("** In Modify T_Mat: Storey Type   \n");  */
       return(T);
    }

    /* general case */

    for(j = 1; j<= no_dof;j++){
        ndof = frp->rigid[rb_no -1].rest_dof[j];
        if(ndof >= 1) { /* dofs not connected by Rigid body are = 1 */
           for(k = 1; k <= no_dof; k++)
               T[j][k] =0.0;
               T[j][j] = 1.0;
        }
    }

    return(T);
}

/* 
 *  ============================================
 *  Set Element Properties for Feap problem        
 *  Input   - p array                               
 *  output  - p array                               
 *  ============================================
 */ 

#ifdef __STDC__
ARRAY *Element_Property(ARRAY *p)
#else
ARRAY *Element_Property(p)
ARRAY *p;
#endif
{
   p = elmlib(p,PROPTY);
   return(p);
}

/* ================================================ */
/* Set Material Properties for Feap problem         */
/* Input  - p array                                 */
/* output  - p array                                */
/* ================================================ */

#ifdef __STDC__
ARRAY *Mate_Property(ARRAY *p)
#else
ARRAY *Mate_Property(p)
ARRAY *p;
#endif
{
int isw = 15;

    p = elmlib(p,isw);
    return(p);
}

/* ================================================== */
/* Assemble Element Property Vector                   */
/* Input  - frame pointer frp, prob case              */
/* ================================================== */

#ifdef __STDC__
ARRAY *Eload_Property(ARRAY *p)
#else
ARRAY *Eload_Property(p)
ARRAY *p;