fe_matrix.c

利用语言编写的有限元分析软件

/* ================================================ */
/* 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;
#endif
{
int  isw = 16;

    p = elmlib(p,isw) ;
    return(p);
}
 
#define PLASTIC 2

/* ======================================================*/
/* function  double   *Transform_Force(frp,F,nn )        */
/*  Calculates  Transform_Force Matrix  Ti   for         */
/*  Relation:        Abi = Ti . Ami                      */
/*  Converts the actions in Ami(at end of the member)    */
/*  to the statically equivalent actions in Abi          */
/*         ( at the working points of the rigid bodies ) */
/* Input - frame pointer frp,                            */
/*      F: Force Vector for the problem which is partially */
/*            modified for node nu 'nn' in transforming  */
/*            node force from this point to wkg pt of RB */
/* ===================================================== */

#ifdef __STDC__
QUANTITY *Transform_Force(EFRAME *frp, QUANTITY *F, int nn )
#else
QUANTITY *Transform_Force(frp,F,nn )
EFRAME   *frp;
QUANTITY *F;
int nn ;
#endif
{
int      i,k,rb_no,rb_elmt_type;
double   **Fj, **Ff, **Tfj;
double   cgx,cgy,cgz ,temp;

   k   = frp->no_dof;       /*for general case*/
   Fj  = MatrixAllocIndirectDouble(k,1);
   Tfj = MatrixAllocIndirectDouble(k,k);

   rb_no = frp->node[nn -1].rb_num;
   rb_elmt_type = 7;
   rb_elmt_type = frp->rigid[rb_no -1].rb_type;

   /* ------------------------------------------------------------------- */
   /* compute coordinates of node wrt ref point- cg(wkg pt) of rigid body */
   /* ------------------------------------------------------------------- */

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

/********
    if(frp->ndm >=3)
       cgz  =   frp->node[nn-1].coord[3] - frp->rigid[rb_no -1].zcg;
    else
       cgz  = 0.0;

**********/
    /* Transformation Matrix */

    Tfj = (double **)Transformation_Matrix(Tfj,(char *)rb_elmt_type,frp->no_dof, cgx,cgy,cgz);
    Tfj = (double **)Modify_T_Matrix(Tfj, frp,frp->no_dof ,rb_no);

    /* Transformation Matrix- store for further use */

    frp->node[nn -1].TrT->uMatrix.daa = (double **) dMatrixTranspose(Tfj,frp->no_dof,frp->no_dof);

    for(i = 1; i<=k;i++)
        Fj[i][1] = F[(nn  - 1) * frp->no_dof +i].value;

    Ff = (double **) dMatrixMult(Tfj, frp->no_dof, frp->no_dof,Fj, frp->no_dof, 1);
    for(i = 1; i<=k;i++)
        F[(nn  - 1) * frp->no_dof+i].value = Ff[i][1];

    MatrixFreeIndirectDouble(Ff, frp->no_dof);
    MatrixFreeIndirectDouble(Fj, frp->no_dof);
    MatrixFreeIndirectDouble(Tfj, frp->no_dof);

    return(F);
}


/* 
 *  ============================================================================ 
 *  Bound_Disp : Checks for displacement at restraint boundary nodes for an elmt             
 *  
 *  Input - frame pointer frp, element no  
 *  ============================================================================ 
 */ 

#ifdef __STDC__
int Bound_Disp(EFRAME *frp, ARRAY *p, int elmt_no)
#else
int Bound_Disp(frp,p, elmt_no)
EFRAME  *frp;
ARRAY   *p;
int     elmt_no;
#endif
{
ELEMENT *el;
double  displ;
int     i,j,k,nen, no_dof,out;

    el  = &frp->element[elmt_no -1];      /* element ptr   */
    p   = Assign_p_Array(frp,elmt_no,p, 10);
    nen = p->nodes_per_elmt;
    no_dof = p->dof_per_node;   /*ed feb 6  MIN(p->dof_per_node,frp->no_dof)*/

    out = NO;         /* default outcome  initialized */

    for(i = 1; i<=nen;i++){
        for(j = 1; j<=no_dof;j++){
            k = el->d_array[(i-1) *no_dof + j];
            if(k<0) {           /* fixed or restrained dofs */
               displ = frp->node[el->node_connect[i]  - 1].disp[j].value; 
               p->nodal_loads[(i-1) * no_dof + j].value = displ;
               if(displ != 0.0) out = YES;
            }
            else
               p->nodal_loads[(i-1) * no_dof + j].value = 0.0;
        }
    }

    return(out);
}

/* =========== */
/* Modify_Load */
/* =========== */

#ifdef __STDC__
QUANTITY *Modify_Load(QUANTITY *b, EFRAME *frp, ARRAY *p, int elmt_no)
#else
QUANTITY *Modify_Load(b, frp, p, elmt_no)
QUANTITY *b;
EFRAME   *frp;
ARRAY    *p;
int      elmt_no;
#endif
{
ELEMENT  *el;
MATRIX   *Ke;
int      *Ld;
int      i,j,k;

    el = &frp->element[elmt_no -1];
    p  = Assign_p_Array(frp,elmt_no,p, LOAD_MATRIX);
    Ke = Element_Matrix(p,STIFF);
    Ld = el->d_array;

    for(i = 1; i<=p->size_of_stiff;i++) {
        k = Ld[i];
        if(k>0) {
           for(j = 1; j<=p->size_of_stiff;j++)
               b[k].value = b[k].value - Ke->uMatrix.daa[i][j] * p->nodal_loads[j].value;
        }
    }

    return(b);
}


/* ============================================== */
/*   Function to Add input loadvector pl          */
/*                      to loadvector fvector     */
/* ============================================== */

#ifdef __STDC__
QUANTITY *Addload_Vector(EFRAME *frp, ARRAY *p, QUANTITY *fvector,  QUANTITY *pl, int elmt_no)
#else
QUANTITY *Addload_Vector(frp,p,fvector, pl, elmt_no)
EFRAME   *frp;
ARRAY    *p;
QUANTITY *fvector, *pl;
int      elmt_no;
#endif
{
int      j,k;
int      i,iid,ii;
ELEMENT  *el;
int      *ldofg;
char     *el_type;
int      elmt_attr_no;

    el           = &frp->element[elmt_no -1];
    elmt_attr_no = el->elmt_attr_no;
    el_type      = frp->eattr[elmt_attr_no-1].elmt_type;
    ldofg        = frp->eattr[elmt_no-1].map_ldof_to_gdof;

    for(i=1;i<=p->nodes_per_elmt;i++){
        ii  = el->node_connect[i];
        iid = (ii-1) * frp->no_dof;
        for(j = 1; j<= p->dof_per_node; j++){
            k = ldofg[j-1];
            if(k > 0)
              fvector[iid + k].value += pl[(i-1)*p->dof_per_node+ j].value;
        }
    }

    return(fvector);
}

/* =============================== */
/* Function to Assemble Nodal Load */
/* =============================== */

#ifdef __STDC__
QUANTITY *Assemble_Nodal_Load(EFRAME *frp, QUANTITY *fv)
#else
QUANTITY *Assemble_Nodal_Load(frp, fv)
EFRAME   *frp;
QUANTITY *fv;
#endif
{
NODE               *np;
NODE_LOADS        *nlp;
int i,node, dof, nload;

 /* Initialize load vector */

    for(dof=1; dof<=frp->no_eq; dof++)
        fv[dof].value = 0.;

    for(nload = 1; nload <= frp->no_node_loads; nload++) {
        nlp  = &frp->nforces[nload-1];
        node = nlp->node_f;
        for(i =1; i <= frp->no_dof; i++)
            fv[(node -1) *  frp->no_dof + i ].value = nlp->fn[i].value;
    }

    return (fv);
}

/* ======================================= */
/*   Function to Assemble  Gravity  Load   */
/* ======================================= */

#ifdef __STDC__
QUANTITY *Assemble_Gravity_Load(EFRAME *frp, QUANTITY *fv)
#else
QUANTITY *Assemble_Gravity_Load(frp, fv)
EFRAME   *frp;
QUANTITY *fv;
#endif
{
NODE               *np;
NODE_LOADS        *nlp;
int i,node, dof, nload;

   /* FILL IN DETAILS LATER */

   return (fv);
}

/* =========================================== */
/*   Function to Apply the Boundary conditions */
/*   to the global stiffness matrix            */
/* =========================================== */

#ifdef __STDC__
double **Boundary_Conditions(EFRAME *frp)
#else
double **Boundary_Conditions(frp)
EFRAME *frp;
#endif
{
NODE                *np;
NODE_LOADS         *nlp;
int  i,node, dof, nload;

   for(node=1;node<=frp->no_nodes;node++) {
       np   = &frp->node[node-1];
       for(i =1; i<= frp->no_dof; i++) {
       if(np->bound_id[i] >0) {
           for(dof=1;dof<=TDOF;dof++) {
               K[dof][(node-1)*frp->no_dof +i] = 0.;
               K[(node-1)*frp->no_dof + 1][dof] = 0.;
           }

           K[(node-1)*frp->no_dof + 1][(node-1)*frp->no_dof+ 1] = 0.;
           F[(node-1)*frp->no_dof + 1]   = 0.;
       }
       }
   }

   return (K);
}

/* ==================================== */
/*   Function to compute nodal forces   */
/* ==================================== */

#ifdef __STDC__
QUANTITY *Nodal_Forces(EFRAME *frp)
#else
QUANTITY *Nodal_Forces(frp)
EFRAME   *frp;
#endif
{
int i,j;

    /* Initialize nodal force vector : Compute [K]*{d} = {Fn} */

     for(i=1; i<=TDOF; i++) Fn[i].value = 0.;
         for(i=1;i<=TDOF;i++)
         for(j=1;j<=TDOF;j++)
             Fn[i].value += K[i][j]*d[j];

     return (Fn);
}

/* =======================================  */
/*   Function to Assemble  Centrifugal  Load*/
/* =======================================  */

#ifdef __STDC__
QUANTITY *Assemble_Ctrfgl_Load(EFRAME *frp, QUANTITY *fv)
#else
QUANTITY *Assemble_Ctrfgl_Load(frp,fv)
EFRAME   *frp;
QUANTITY *fv;
#endif
{
NODE               *np;
NODE_LOADS        *nlp;
int i,node, dof, nload;

 /* to be changed */

   for(dof=1; dof<=frp->no_eq; dof++)
      fv[dof].value   = 0.;

   for(nload = 1; nload <= frp->no_node_loads; nload++) {

     nlp  = &frp->nforces[nload-1];
     node = nlp->node_f;
         for(i =1; i <= frp->no_dof; i++)
             fv[(node -1) *  frp->no_dof + i ].value = nlp->fn[i].value;

   }

   return (fv);
}



/* Functions for design rule checking */

#ifdef __STDC__
MATRIX *Get_Coord(MATRIX *m)
#else
MATRIX *Get_Coord(m)
MATRIX *m;
#endif
{
MATRIX *coord;
int    nodeno;
int        ii;

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

       nodeno = (int) m->uMatrix.daa[0][0];
       coord  = MatrixAllocIndirect("Node Coord", DOUBLE_ARRAY, 1, frame->no_dimen);

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