fe_matrix.c

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

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *          fe_matrix.c : Functions to solve (non)linear FE solution procedures
 *                                                                     
 *  Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
 *  Institute for Systems Research,                                           
 *  University of Maryland, College Park, MD 20742                                   
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software for any on any computer system
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *                                                                    
 *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
 *  ============================================================================= 
 */

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#ifdef  __STDC__
#include <stdarg.h>
#else
#include <varargs.h>
#endif

#include "defs.h"
#include "units.h"
#include "matrix.h"
#include "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"
#include "vector.h"
#include "elmt.h"

extern ARRAY     *array;
extern EFRAME    *frame;

static QUANTITY     *Fn;

/*
#define DEBUG
*/

/* 
 *  ------------------------------ 
 *  Form Stiffness and Mass Matrix
 *  ------------------------------ 
 */ 

MATRIX *Form_Stiffness() {
MATRIX      *stiff;
ELEMENT        *el;
MATRIX         *Ke; 
int            *Ld;
int        elmt_no;
int            i,j;
int        iMinRow;
int      iColumnNo;
int   *ipColHeight;

   /* [a] : Compute skyline profile for global stiffness matrix */

   ipColHeight = (int *)iVectorAlloc( TNEQ );

   for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
       Ld    = Destination_Array(frame, elmt_no);
       j = 1;
       while( Ld[j]<0 ) j++;
       iMinRow = Ld[j];
       for( i=j+1 ; i<=Ld[0] ; i++ ) {
            if( Ld[i]>0 )
                iMinRow   = MIN(iMinRow, Ld[i]);
       }

       for( i=j ; i<=Ld[0] ; i++ ) {
            if( Ld[i]>0 )  iColumnNo = Ld[i];
            else           iColumnNo = 1;
            if((1+iColumnNo-iMinRow) > ipColHeight[iColumnNo-1])
                ipColHeight[iColumnNo-1] = 1+iColumnNo-iMinRow;
       }
       free((char *) Ld);
   }

   /* [b] : Allocate memory for global stiffness matrix */

    stiff = MatrixAllocSkyline("stiff", DOUBLE_ARRAY, TNEQ, TNEQ, ipColHeight);
    free((char *)ipColHeight);

   /* [c] : Add element level stiffness matrices into global stiffness matrix */

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

       array = Assign_p_Array(frame, elmt_no, array, STIFF);
       array = Assign_p_Array(frame, elmt_no, array, PROPTY);

       array = Element_Property(array);  /* passing elmt info in             */
                                         /* elmt liberary to working array   */
       Ke = Element_Matrix(array, STIFF);/* set elment stiffness/mass matrix */

       if(Rigidbody_conn(frame,array) == YES)
          Ke = Transform_Stiff_Matrix(frame,array,Ke);

       /* transfer Destination array from frame to Ld */

       Ld    = Destination_Array(frame, elmt_no); 
        
       stiff = Assemble_Global(stiff, frame, Ld, Ke);

       free((char *) Ld);
    } 

    return(stiff);
}

/*
 *  ========================================================
 *  Form Equivalent Nodal Load due to body force      
 *  inital stress, initial strain and surface loadings
 *  ========================================================
 */

MATRIX *Form_Equiv_Nodal_Load() {
MATRIX  *equiv_nodal_load;
MATRIX        *nodal_load;
ELEMENT               *el;
MATRIX                *Fe; 
int                   *Ld;
int               elmt_no;
int          dof_per_elmt;
int                   i,j;

   /* [a] : Allocate memory for equivalent nodal load vector */

    dof_per_elmt     = frame->no_nodes_per_elmt*frame->no_dof;
    equiv_nodal_load = MatrixAllocIndirect("equiv_nodal_load", DOUBLE_ARRAY, TNEQ, 1);

   /* [b] : Add element level nodal loads into global load matrix */

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

        array = Assign_p_Array(frame, elmt_no, array, EQUIV_NODAL_LOAD);
        Fe    = Element_Equiv(array, EQUIV_NODAL_LOAD);
        Ld    = Destination_Array(frame, elmt_no); 

        equiv_nodal_load = Assemble_Global_Load(equiv_nodal_load, frame, Ld, Fe);
    }

    return(equiv_nodal_load);
}

/* 
 *  =================================================
 *  Form Global Mass Matrix
 *  
 *  Input  : frame, p, type = CONSISTENT or LUMPED.
 *  Output : MATRIX    *mass
 *  =================================================
 */ 

#ifdef __STDC__
MATRIX *Form_Mass(MATRIX *m)
#else
MATRIX *Form_Mass(m)
MATRIX *m;
#endif
{
MATRIX       *mass;
ELEMENT        *el;
MATRIX         *Me; 
int            *Ld;
int        elmt_no;
int       rigid_no;
int      mass_type;
int       size,i,j;
int        iMinRow;
int      iColumnNo;
int   *ipColHeight;

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

   /* [a] : Setup type of mass[][] to be computed */ 

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

   /* [b] : Compute skyline profile for mass matrix */

   ipColHeight = (int *)iVectorAlloc( TNEQ );

   for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
       Ld    = Destination_Array(frame, elmt_no);
       j = 1;
       while( Ld[j]<0 ) j++;

       iMinRow = Ld[j];
       for( i=j+1 ; i<=Ld[0] ; i++ ) {
            if( Ld[i]>0 )
                iMinRow   = MIN(iMinRow, Ld[i]);
       }

       for( i=j ; i<=Ld[0] ; i++ ) {
            if( Ld[i]>0 )  iColumnNo = Ld[i];
            else           iColumnNo = 1;
            if((1+iColumnNo-iMinRow) > ipColHeight[iColumnNo-1])
                ipColHeight[iColumnNo-1] = 1+iColumnNo-iMinRow;
       }
       free((char *) Ld);
   }

   /* [c] : Allocate memory for skyline mass matrix */

   mass = MatrixAllocSkyline("mass", DOUBLE_ARRAY, TNEQ, TNEQ, ipColHeight);
   free((char *)ipColHeight);

   /* [d] : Form global mass matrix from element mass matrices */

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

      array = Assign_p_Array(frame, elmt_no, array, MASS_MATRIX);
      array = Assign_p_Array(frame, elmt_no, array, PROPTY);
      array->type = mass_type;

      array = Element_Property(array);

      Me = Element_Matrix(array, MASS_MATRIX);

      /* If element connected to rigid body transform  */
      /* mass matrix to working points of Rigid Body   */

      if(Rigidbody_conn(frame,array) == YES)
         Me = Transform_Stiff_Matrix(frame,array,Me);

      /* Me : Local element stiffness/mass matrix, size_of_stiffxsize_of_stiff */
      /* mass : Global stiffness/mass matrix, TNEQxTNEQ    */

      Ld   = Destination_Array(frame, elmt_no);
      mass = Assemble_Global(mass, frame, Ld, Me);

      free((char *) Ld);
   }

   if(frame->no_rigid == 0) {
      return (mass);
   }

   /*
    *  -----------------------------------------------------------------
    *  [e] : Compute Mass Matrix for Rigid Body Components
    * 
    *        Calculate [m] for rigid bodies; add to global mass matrix.
    * 
    *        Transform mass matrix from mass center 'c' to working pt 'p'
    *        by : Mp[][] = Tpc[][].Mc[][].Ttpc[][].
    *  -----------------------------------------------------------------
    */

    /*  ----- fix this section of code later -----------------------

    for(rigid_no = 1; rigid_no <= frame->no_rigid; rigid_no++) {

        array       = Assign_p_Array_for_Rigid_Body(frame, rigid_no, array, MASS_MATRIX);
        array->type = mass_type;                                                          * LUMPED or CONSISTENT * 

        Me = Element_Matrix(array,MASS_MATRIX);
        Me = Transform_Rigid_Body_Mass_Matrix(frame, array, Me);

        Ld   = Destination_Array_for_Rigid_Body(frame, array, rigid_no, NO);
        mass = Assemble_Global(mass, frame, Ld, Me);
    }

     ---- fix this block of code later --------------------------*/

    return(mass);
}

/* 
 *  =======================================
 *  Form External and Internal Load Vectors
 *  =======================================
 */ 

MATRIX *Form_External_Load() {
MATRIX                     *load;
NODE_LOADS              *nforces;
DIMENSIONS            *units_buf;
int node_no, i, j, k,jj, counter;
int node_no_1,            length;
int                 UNITS_SWITCH;

#ifdef DEBUG
       printf("*** Enter Form_External_Load() TNEQ = %d\n", TNEQ);
#endif

   /* [a] : Allocate memory for external load vector */

   UNITS_SWITCH = CheckUnits();
   load = MatrixAllocIndirect("load", DOUBLE_ARRAY, TNEQ, 1);    

   if( UNITS_SWITCH == ON ) {
       for (i = 1; i <= TNEQ; i++)
            ZeroUnits(&(load->spRowUnits[i-1]));
   }

   /* [b] : Transfer units into units buffer */

   k = 1;
   nforces = &frame->nforces[k-1];
   node_no_1 = nforces->node_f;

   if( UNITS_SWITCH == ON ) {
       for(i = 1; i <= frame->no_nodes; i++) {
           for(j = 1; j <= frame->no_dof; j++) {
               jj = frame->node[i-1].bound_id[j-1];
               if(jj > 0)
                  UnitsCopy( &(load->spRowUnits[jj-1]), nforces->fn[j-1].dimen );
           }
       }
   }

   for(k = 1; k <= frame->no_node_loads; k++) {
       nforces = &frame->nforces[k-1];
       node_no = nforces->node_f;
       for(j = 1; j <= frame->no_dof; j++) {
           jj = frame->node[node_no-1].bound_id[j-1];
           if(jj > 0)
              load->uMatrix.daa[(int)(jj-1)][0] += nforces->fn[j-1].value;
       }
   }

   if( UNITS_SWITCH == ON ) {
       ZeroUnits(&(load->spColUnits[0]));
       load->spColUnits[0].units_type = nforces->fn[0].dimen->units_type;
   }

#ifdef DEBUG
       MatrixPrintIndirectDouble(load);
       printf("*** Leave Form_External_Load()\n");
#endif

    return(load);
}

/* 
 *  ==============================================================
 *  Form Internal Load Vector
 * 
 *  Usage : load = Form_Internal_Load(displ) 
 *       or load = Form_Internal_Load(displ, displ_incr)         
 *  
 *  Boundary Info : INPUT -- bound_id = 0, no restraint,                
 *                           bound_id > 0, restrained, displ specified  
 *                           bound_id is changed based on INPUT bound_id
 *               : OUTPUT -- bound_id < 0, restrained displ specified    
 *                           bound_id > 0, no restraint                  
 * 
 *  ==============================================================
 */ 
 
#ifdef  __STDC__
MATRIX *Form_Internal_Load( MATRIX *displ, ... )
#else
MATRIX *Form_Internal_Load(va_alist)
va_dcl
#endif
{
#ifndef __STDC__
MATRIX            *displ;
#endif

va_list          arg_ptr;
MATRIX       *displ_incr;
MATRIX             *load;
ELEMENT              *ep;
int     node_no, elmt_no;
int     i, j, k, jj, *Ld; 
int         elmt_attr_no;
int     length1, length2;
int         UNITS_SWITCH;
int itemp;

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

   UNITS_SWITCH = CheckUnits();

   /* [a] : read input */

#ifdef __STDC__
   va_start(arg_ptr, displ);
#else
   va_start(arg_ptr);
   displ = va_arg(arg_ptr, MATRIX *);
#endif

   if(array->material_name != NULL && strcmp(array->material_name, "ELASTIC")) {