fe_matrix.c

矩阵运算〉〉 〉〉 〉〉 〉〉 。

   1: /*
   2:  *  ============================================================================= 
   3:  *  ALADDIN Version 1.0 :
   4:  *          fe_matrix.c : Functions to solve (non)linear FE solution procedures
   5:  *                                                                     
   6:  *  Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
   7:  *  Institute for Systems Research,                                           
   8:  *  University of Maryland, College Park, MD 20742                                   
   9:  *                                                                     
  10:  *  This software is provided "as is" without express or implied warranty.
  11:  *  Permission is granted to use this software for any on any computer system
  12:  *  and to redistribute it freely, subject to the following restrictions:
  13:  * 
  14:  *  1. The authors are not responsible for the consequences of use of
  15:  *     this software, even if they arise from defects in the software.
  16:  *  2. The origin of this software must not be misrepresented, either
  17:  *     by explicit claim or by omission.
  18:  *  3. Altered versions must be plainly marked as such, and must not
  19:  *     be misrepresented as being the original software.
  20:  *  4. This notice is to remain intact.
  21:  *                                                                    
  22:  *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
  23:  *  ============================================================================= 
  24:  */
  25: 
  26: #include <stdio.h>
  27: #include <stdlib.h>
  28: #include <ctype.h>
  29: #ifdef  __STDC__
  30: #include <stdarg.h>
  31: #else
  32: #include <varargs.h>
  33: #endif
  34: 
  35: #include "defs.h"
  36: #include "units.h"
  37: #include "matrix.h"
  38: #include "fe_database.h"
  39: #include "symbol.h"
  40: #include "fe_functions.h"
  41: #include "vector.h"
  42: #include "elmt.h"
  43: 
  44: extern ARRAY     *array;
  45: extern EFRAME    *frame;
  46: 
  47: static QUANTITY     *Fn;
  48: 
  49: /*
  50: #define DEBUG
  51: */
  52: 
  53: /* 
  54:  *  ------------------------------ 
  55:  *  Form Stiffness and Mass Matrix
  56:  *  ------------------------------ 
  57:  */ 
  58: 
  59: MATRIX *Form_Stiffness() {
  60: MATRIX      *stiff;
  61: ELEMENT        *el;
  62: MATRIX         *Ke; 
  63: int            *Ld;
  64: int        elmt_no;
  65: int            i,j;
  66: int        iMinRow;
  67: int      iColumnNo;
  68: int   *ipColHeight;
  69: 
  70:    /* [a] : Compute skyline profile for global stiffness matrix */
  71: 
  72:    ipColHeight = (int *)iVectorAlloc( TNEQ );
  73: 
  74:    for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
  75:        Ld    = Destination_Array(frame, elmt_no);
  76:        j = 1;
  77:        while( Ld[j]<0 ) j++;
  78:        iMinRow = Ld[j];
  79:        for( i=j+1 ; i<=Ld[0] ; i++ ) {
  80:             if( Ld[i]>0 )
  81:                 iMinRow   = MIN(iMinRow, Ld[i]);
  82:        }
  83: 
  84:        for( i=j ; i<=Ld[0] ; i++ ) {
  85:             if( Ld[i]>0 )  iColumnNo = Ld[i];
  86:             else           iColumnNo = 1;
  87:             if((1+iColumnNo-iMinRow) > ipColHeight[iColumnNo-1])
  88:                 ipColHeight[iColumnNo-1] = 1+iColumnNo-iMinRow;
  89:        }
  90:        free((char *) Ld);
  91:    }
  92: 
  93:    /* [b] : Allocate memory for global stiffness matrix */
  94: 
  95:     stiff = MatrixAllocSkyline("stiff", DOUBLE_ARRAY, TNEQ, TNEQ, ipColHeight);
  96:     free((char *)ipColHeight);
  97: 
  98:    /* [c] : Add element level stiffness matrices into global stiffness matrix */
  99: 
 100:    for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
 101: 
 102:        array = Assign_p_Array(frame, elmt_no, array, STIFF);
 103:        array = Assign_p_Array(frame, elmt_no, array, PROPTY);
 104: 
 105:        array = Element_Property(array);  /* passing elmt info in             */
 106:                                          /* elmt liberary to working array   */
 107:        Ke = Element_Matrix(array, STIFF);/* set elment stiffness/mass matrix */
 108: 
 109:        if(Rigidbody_conn(frame,array) == YES)
 110:           Ke = Transform_Stiff_Matrix(frame,array,Ke);
 111: 
 112:        /* transfer Destination array from frame to Ld */
 113: 
 114:        Ld    = Destination_Array(frame, elmt_no); 
 115:         
 116:        stiff = Assemble_Global(stiff, frame, Ld, Ke);
 117: 
 118:        free((char *) Ld);
 119:     } 
 120: 
 121:     return(stiff);
 122: }
 123: 
 124: /*
 125:  *  ========================================================
 126:  *  Form Equivalent Nodal Load due to body force      
 127:  *  inital stress, initial strain and surface loadings
 128:  *  ========================================================
 129:  */
 130: 
 131: MATRIX *Form_Equiv_Nodal_Load() {
 132: MATRIX  *equiv_nodal_load;
 133: MATRIX        *nodal_load;
 134: ELEMENT               *el;
 135: MATRIX                *Fe; 
 136: int                   *Ld;
 137: int               elmt_no;
 138: int          dof_per_elmt;
 139: int                   i,j;
 140: 
 141:    /* [a] : Allocate memory for equivalent nodal load vector */
 142: 
 143:     dof_per_elmt     = frame->no_nodes_per_elmt*frame->no_dof;
 144:     equiv_nodal_load = MatrixAllocIndirect("equiv_nodal_load", DOUBLE_ARRAY, TNEQ, 1);
 145: 
 146:    /* [b] : Add element level nodal loads into global load matrix */
 147: 
 148:     for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
 149: 
 150:         array = Assign_p_Array(frame, elmt_no, array, EQUIV_NODAL_LOAD);
 151:         Fe    = Element_Equiv(array, EQUIV_NODAL_LOAD);
 152:         Ld    = Destination_Array(frame, elmt_no); 
 153: 
 154:         equiv_nodal_load = Assemble_Global_Load(equiv_nodal_load, frame, Ld, Fe);
 155:     }
 156: 
 157:     return(equiv_nodal_load);
 158: }
 159: 
 160: /* 
 161:  *  =================================================
 162:  *  Form Global Mass Matrix
 163:  *  
 164:  *  Input  : frame, p, type = CONSISTENT or LUMPED.
 165:  *  Output : MATRIX    *mass
 166:  *  =================================================
 167:  */ 
 168: 
 169: #ifdef __STDC__
 170: MATRIX *Form_Mass(MATRIX *m)
 171: #else
 172: MATRIX *Form_Mass(m)
 173: MATRIX *m;
 174: #endif
 175: {
 176: MATRIX       *mass;
 177: ELEMENT        *el;
 178: MATRIX         *Me; 
 179: int            *Ld;
 180: int        elmt_no;
 181: int       rigid_no;
 182: int      mass_type;
 183: int       size,i,j;
 184: int        iMinRow;
 185: int      iColumnNo;
 186: int   *ipColHeight;
 187: 
 188: #ifdef DEBUG
 189:        printf("*** Enter Form_Mass()\n");
 190: #endif
 191: 
 192:    /* [a] : Setup type of mass[][] to be computed */ 
 193: 
 194:    mass_type = (int) m->uMatrix.daa[0][0];
 195: 
 196:    /* [b] : Compute skyline profile for mass matrix */
 197: 
 198:    ipColHeight = (int *)iVectorAlloc( TNEQ );
 199: 
 200:    for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
 201:        Ld    = Destination_Array(frame, elmt_no);
 202:        j = 1;
 203:        while( Ld[j]<0 ) j++;
 204: 
 205:        iMinRow = Ld[j];
 206:        for( i=j+1 ; i<=Ld[0] ; i++ ) {
 207:             if( Ld[i]>0 )
 208:                 iMinRow   = MIN(iMinRow, Ld[i]);
 209:        }
 210: 
 211:        for( i=j ; i<=Ld[0] ; i++ ) {
 212:             if( Ld[i]>0 )  iColumnNo = Ld[i];
 213:             else           iColumnNo = 1;
 214:             if((1+iColumnNo-iMinRow) > ipColHeight[iColumnNo-1])
 215:                 ipColHeight[iColumnNo-1] = 1+iColumnNo-iMinRow;
 216:        }
 217:        free((char *) Ld);
 218:    }
 219: 
 220:    /* [c] : Allocate memory for skyline mass matrix */
 221: 
 222:    mass = MatrixAllocSkyline("mass", DOUBLE_ARRAY, TNEQ, TNEQ, ipColHeight);
 223:    free((char *)ipColHeight);
 224: 
 225:    /* [d] : Form global mass matrix from element mass matrices */
 226: 
 227:    for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
 228: 
 229:       array = Assign_p_Array(frame, elmt_no, array, MASS_MATRIX);
 230:       array = Assign_p_Array(frame, elmt_no, array, PROPTY);
 231:       array->type = mass_type;
 232: 
 233:       array = Element_Property(array);
 234: 
 235:       Me = Element_Matrix(array, MASS_MATRIX);
 236: 
 237:       /* If element connected to rigid body transform  */
 238:       /* mass matrix to working points of Rigid Body   */
 239: 
 240:       if(Rigidbody_conn(frame,array) == YES)
 241:          Me = Transform_Stiff_Matrix(frame,array,Me);
 242: 
 243:       /* Me : Local element stiffness/mass matrix, size_of_stiffxsize_of_stiff */
 244:       /* mass : Global stiffness/mass matrix, TNEQxTNEQ    */
 245: 
 246:       Ld   = Destination_Array(frame, elmt_no);
 247:       mass = Assemble_Global(mass, frame, Ld, Me);
 248: 
 249:       free((char *) Ld);
 250:    }
 251: 
 252:    if(frame->no_rigid == 0) {
 253:       return (mass);
 254:    }
 255: 
 256:    /*
 257:     *  -----------------------------------------------------------------
 258:     *  [e] : Compute Mass Matrix for Rigid Body Components
 259:     * 
 260:     *        Calculate [m] for rigid bodies; add to global mass matrix.
 261:     * 
 262:     *        Transform mass matrix from mass center 'c' to working pt 'p'
 263:     *        by : Mp[][] = Tpc[][].Mc[][].Ttpc[][].
 264:     *  -----------------------------------------------------------------
 265:     */
 266: 
 267:     /*  ----- fix this section of code later -----------------------
 268: 
 269:     for(rigid_no = 1; rigid_no <= frame->no_rigid; rigid_no++) {
 270: 
 271:         array       = Assign_p_Array_for_Rigid_Body(frame, rigid_no, array, MASS_MATRIX);
 272:         array->type = mass_type;                                                          * LUMPED or CONSISTENT * 
 273: 
 274:         Me = Element_Matrix(array,MASS_MATRIX);
 275:         Me = Transform_Rigid_Body_Mass_Matrix(frame, array, Me);
 276: 
 277:         Ld   = Destination_Array_for_Rigid_Body(frame, array, rigid_no, NO);
 278:         mass = Assemble_Global(mass, frame, Ld, Me);
 279:     }
 280: 
 281:      ---- fix this block of code later --------------------------*/
 282: 
 283:     return(mass);
 284: }
 285: 
 286: /* 
 287:  *  =======================================
 288:  *  Form External and Internal Load Vectors
 289:  *  =======================================
 290:  */ 
 291: 
 292: MATRIX *Form_External_Load() {
 293: MATRIX                     *load;
 294: NODE_LOADS              *nforces;
 295: DIMENSIONS            *units_buf;
 296: int node_no, i, j, k,jj, counter;
 297: int node_no_1,            length;
 298: int                 UNITS_SWITCH;
 299: 
 300: #ifdef DEBUG
 301:        printf("*** Enter Form_External_Load() TNEQ = %d\n", TNEQ);
 302: #endif
 303: 
 304:    /* [a] : Allocate memory for external load vector */
 305: 
 306:    UNITS_SWITCH = CheckUnits();
 307:    load = MatrixAllocIndirect("load", DOUBLE_ARRAY, TNEQ, 1);    
 308: 
 309:    if( UNITS_SWITCH == ON ) {
 310:        for (i = 1; i <= TNEQ; i++)
 311:             ZeroUnits(&(load->spRowUnits[i-1]));
 312:    }
 313: 
 314:    /* [b] : Transfer units into units buffer */
 315: 
 316:    k = 1;
 317:    nforces = &frame->nforces[k-1];
 318:    node_no_1 = nforces->node_f;
 319: 
 320:    if( UNITS_SWITCH == ON ) {
 321:        for(i = 1; i <= frame->no_nodes; i++) {
 322:            for(j = 1; j <= frame->no_dof; j++) {
 323:                jj = frame->node[i-1].bound_id[j-1];
 324:                if(jj > 0)
 325:                   UnitsCopy( &(load->spRowUnits[jj-1]), nforces->fn[j-1].dimen );
 326:            }
 327:        }
 328:    }
 329: 
 330:    for(k = 1; k <= frame->no_node_loads; k++) {
 331:        nforces = &frame->nforces[k-1];
 332:        node_no = nforces->node_f;
 333:        for(j = 1; j <= frame->no_dof; j++) {
 334:            jj = frame->node[node_no-1].bound_id[j-1];
 335:            if(jj > 0)
 336:               load->uMatrix.daa[(int)(jj-1)][0] += nforces->fn[j-1].value;
 337:        }
 338:    }
 339: 
 340:    if( UNITS_SWITCH == ON ) {
 341:        ZeroUnits(&(load->spColUnits[0]));
 342:        load->spColUnits[0].units_type = nforces->fn[0].dimen->units_type;
 343:    }
 344: 
 345: #ifdef DEBUG
 346:        MatrixPrintIndirectDouble(load);
 347:        printf("*** Leave Form_External_Load()\n");
 348: #endif
 349: 
 350:     return(load);
 351: }
 352: 
 353: /* 
 354:  *  ==============================================================
 355:  *  Form Internal Load Vector