elmt_fiber2d.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  elmt_fiber2d.c : Linear/Nonlinear 2D Fiber Element
 *                                                                     
 *  Copyright (C) 1995-2000 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 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 software may not be sold or included in commercial software
 *     products without a license. 
 *  5. This notice is to remain intact.
 *                                                                    
 *  --------------------------------------------------------------------
 *  Convention for Nodal Forces
 *             +ve M     -  anticlockwise(RHT Rule)
 *             +ve X,Y,Z -  along +ve axis
 *  Convention for Member End  Forces
 *             +ve M     -  Sagging Moments
 *             +ve SF    -  LHS upwards
 *             +ve AF    -  Tension(LHS outwards)
 *  --------------------------------------------------------------------
 *  Written by: Wane-Jang Lin                                   May 1996
 *  Modified by: Mark Austin                                  March 2000
 *  ====================================================================
 */

#include <math.h>
#include "defs.h"
#include "miscellaneous.h"
#include "units.h"
#include "matrix.h"
#include "fe_database.h"
#include "fe_functions.h"
#include "elmt.h"

/*
#define  DEBUG 
*/

/* 
 *  ============================================================== 
 *  Element FIBER_2D                                            
 *
 *  Input Properties:                                      
 *
 *       p->work_material[0] = E;
 *       p->work_material[1] = G;
 *       p->work_material[2] = fy;
 *       p->work_material[3] = ET;
 *       p->work_material[4] = nu;
 *       p->work_material[5] = density;
 *       p->work_material[6] = fu;
 *
 *       p->work_section[0] = Ixx;
 *       p->work_section[1] = Iyy;
 *       p->work_section[2] = Izz;
 *       p->work_section[3] = Ixy;
 *       p->work_section[4] = Ixz;
 *       p->work_section[5] = Iyz;
 *       p->work_section[6] = weight;
 *       p->work_section[7] = bf;
 *       p->work_section[8] = tf;
 *       p->work_section[9] = depth;                                  
 *       p->work_section[10] = area;
 *       p->work_section[11] = plate_thickness;
 *       p->work_section[12] = J;
 *       p->work_section[13] = rT;
 *       p->work_section[14] = width;
 *       p->work_section[15] = tw;                                  
 * 
 *  Input  :  ARRAY *p  -- pointer to working ARRAY data structure
 *         :  int isw   -- flag for task to be computed.
 *  Output :  ARRAY *p  -- pointer to working ARRAY data structure
 *  ============================================================== 
 */

#ifdef __STDC__
ARRAY *elmt_fiber_2d(ARRAY *p, int isw)
#else
ARRAY *elmt_fiber_2d(p, isw)
ARRAY *p;
int     isw;
#endif
{
static int       no_integ_pt, no_section, no_fiber, no_shear, elmt_no;
static QUANTITY  elmt_length;

int         ii, jj, kk, sec, ifib;
int         UNITS_SWITCH, UnitsType;
double      cs, sn, tn, xx, yy;
double      mbar;
MATRIX      *temp_m1, *temp_m2;
DIMENSIONS  *dp_length, *dp_force, *dp_stress, *dimen;

double      *abscissas, *weights;
double      xi;
double      scale;
MATRIX      *Q, *q;
MATRIX      *F, *K, *Ke;
MATRIX      *L, *Ltrans, *R, *Rtrans;
MATRIX      *kx, *fx;
MATRIX      *bx, *bxtrans;
FIBER_ATTR  *fiber;

    UNITS_SWITCH = CheckUnits();
    UnitsType    = CheckUnitsType();

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

    switch( isw ) {
       case PROPTY :
            no_fiber    = p->fiber_ptr->no_fiber;
            no_shear    = p->fiber_ptr->no_shear;
            no_integ_pt = p->integ_ptr->integ_pts;
            no_section  = no_integ_pt + 2;
            elmt_no     = p->elmt_no;
            cs = p->coord[0][1].value - p->coord[0][0].value;
            sn = p->coord[1][1].value - p->coord[1][0].value;
            elmt_length.value = sqrt( cs*cs + sn*sn );

            if( UNITS_SWITCH == ON ) {
               if( UnitsType == SI )
                  dp_length = DefaultUnits("m");
               else if( UnitsType == US )
                  dp_length = DefaultUnits("in");

               elmt_length.dimen = (DIMENSIONS *)MyCalloc(1, sizeof(DIMENSIONS));
               UnitsCopy( elmt_length.dimen, dp_length );
               free((char *) dp_length->units_name);
               free((char *) dp_length);
            }

            break;  /* end of case PROPTY */
       case STIFF :
            if( UNITS_SWITCH == ON )
                SetUnitsOff();

            fiber     = p->fiber_ptr->fiber;
            abscissas = (double *)MyCalloc( no_section, sizeof(double) );
            weights   = (double *)MyCalloc( no_section, sizeof(double) );
            Gauss_Lobatto( abscissas, weights, no_integ_pt );

            /* Allocate working arrays for element flexibility */

            kx = MatrixAllocIndirect( "kx", DOUBLE_ARRAY, 3, 3 );
            bx = MatrixAllocIndirect( "bx", DOUBLE_ARRAY, 3, 3 );

            /* Compute element flexibility matrix */

            for( sec=0 ; sec < no_section; ++sec ) {
                xi    = elmt_length.value/2.0*(abscissas[sec]+1);
                scale = weights[sec]*elmt_length.value/2.0;

                Force_Interpolation_Matrix_2d( bx, xi, elmt_length );

                Section_Tangent_Stiffness_2d( kx, fiber, no_fiber, no_shear, sec, elmt_no );

                fx      = MatrixInverse( kx );
                bxtrans = MatrixTranspose( bx );
                temp_m1 = MatrixMult( bxtrans, fx );
                temp_m2 = MatrixMult( temp_m1, bx );

                MatrixFree( fx );
                MatrixFree( temp_m1 );
                MatrixFree( bxtrans );

                if( sec == 0 )
                   F = MatrixScale( temp_m2, scale );
                else {
                   temp_m1 = MatrixScale( temp_m2, scale );
                   MatrixAddReplace( F, temp_m1 );
                   MatrixFree( temp_m1 );
                }

                MatrixFree( temp_m2 );
            }

            /* Invert element flexibility to get stiffness matrix */

            K = MatrixInverse( F );

            free((char *) abscissas);
            free((char *) weights);
            MatrixFree( kx );
            MatrixFree( bx );
            MatrixFree( F );

            /* Rigid body rotation and transform local coordinate to global */

            L       = Element_Transformation_2d( p->coord, elmt_length );
            Ltrans  = MatrixTranspose( L );
            temp_m1 = MatrixMult( Ltrans, K );
            Ke      = MatrixMult( temp_m1, L );

            /* Copy stiffness matrix to p array */

            for(ii = 1; ii <= p->stiff->iNoRows; ii++)
            for(jj = 1; jj <= p->stiff->iNoColumns; jj++)
                p->stiff->uMatrix.daa[ii-1][jj-1] = Ke->uMatrix.daa[ii-1][jj-1];

            MatrixFree( temp_m1 );
            MatrixFree( K );
            MatrixFree( L );
            MatrixFree( Ltrans );
            MatrixFree( Ke );

            /* Assign units to p array stiffness */

            if( UNITS_SWITCH == ON ) {
                SetUnitsOn();
                switch( UnitsType ) {
                   case SI:
                   case SI_US:
                        dp_force  = DefaultUnits("N");
                        dp_length = DefaultUnits("m");
                        break;
                   case US:
                        dp_force  = DefaultUnits("lbf");
                        dp_length = DefaultUnits("in");
                        break;
                }

                ZeroUnits( &(p->stiff->spRowUnits[0]) );
                ZeroUnits( &(p->stiff->spRowUnits[1]) );
                UnitsCopy( &(p->stiff->spRowUnits[2]), dp_length );
                UnitsCopy( &(p->stiff->spRowUnits[3]),
                           &(p->stiff->spRowUnits[0]) );
                UnitsCopy( &(p->stiff->spRowUnits[4]),
                           &(p->stiff->spRowUnits[1]) );
                UnitsCopy( &(p->stiff->spRowUnits[5]),
                           &(p->stiff->spRowUnits[2]) );

                UnitsDivRep( &(p->stiff->spColUnits[0]),
                             dp_force, dp_length, NO );
                UnitsCopy( &(p->stiff->spColUnits[1]),
                           &(p->stiff->spColUnits[0]) );
                UnitsCopy( &(p->stiff->spColUnits[2]), dp_force );
                UnitsCopy( &(p->stiff->spColUnits[3]),
                           &(p->stiff->spColUnits[0]) );
                UnitsCopy( &(p->stiff->spColUnits[4]),
                           &(p->stiff->spColUnits[1]) );
                UnitsCopy( &(p->stiff->spColUnits[5]),
                           &(p->stiff->spColUnits[2]) );

                free((char *) dp_force->units_name);
                free((char *) dp_length->units_name);
                free((char *) dp_force);
                free((char *) dp_length);
            } /* end of units on/off for case STIFF */

            break;   /* end of case STIFF */

       case LOAD_MATRIX:  /* compute internal load vector       */
       case STRESS:       /* compute and print element stresses */

            if( UNITS_SWITCH == ON )
               SetUnitsOff();

            /* --------------------------------------------------------- */
            /* Two cases to deal with:                                   */
            /*                                                           */
            /* p->elmt_state == 0 : ElmtStateDet() has not been called.  */
            /*    Therefore, element stresses and loads must be computed */
            /*    from scratch.                                          */
            /* p->elmt_state == 1 : ElmtStateDet() has been called.      */
            /*    Therefore, element stresses and loads are saveed and   */
            /*    only need to be retrieved from Aladdin database.       */
            /* --------------------------------------------------------- */

            if( p->elmt_state == 0 ) { 

            fiber     = p->fiber_ptr->fiber;
            abscissas = (double *)MyCalloc( no_section, sizeof(double) );
            weights   = (double *)MyCalloc( no_section, sizeof(double) );
            Gauss_Lobatto( abscissas, weights, no_integ_pt );

            /* Allocate memory for element force/stiffness */

            kx = MatrixAllocIndirect( "kx", DOUBLE_ARRAY, 3, 3 );
            bx = MatrixAllocIndirect( "bx", DOUBLE_ARRAY, 3, 3 );

            /* Construct element flexibility matrix by looping over sections */

            for( sec=0 ; sec < no_section; ++sec ) {
                xi    = elmt_length.value/2.0*(abscissas[sec]+1);
                scale = weights[sec]*elmt_length.value/2.0;

                Force_Interpolation_Matrix_2d( bx, xi, elmt_length );
                Section_Tangent_Stiffness_2d( kx, fiber, no_fiber, no_shear, sec, elmt_no );

                fx = MatrixInverse( kx );

                bxtrans = MatrixTranspose( bx );
                temp_m1 = MatrixMult( bxtrans, fx );
                temp_m2 = MatrixMult( temp_m1, bx );

                MatrixFree( fx );
                MatrixFree( temp_m1 );
                MatrixFree( bxtrans );

                if( sec == 0 )
                   F = MatrixScale( temp_m2, scale );
                else {
                   temp_m1 = MatrixScale( temp_m2, scale );
                   MatrixAddReplace( F, temp_m1 );
                   MatrixFree( temp_m1 );
                }
                MatrixFree( temp_m2 );
            }

            /* Invert element flexibility to get element-level stiffness */

            K = MatrixInverse( F );

            /* Clean up ..... */

            free((char *) abscissas);
            free((char *) weights);
            MatrixFree( kx );
            MatrixFree( bx );
            MatrixFree( F );

            /* --------------------------------------------------------- */
            /* Put displacement matrix into one column format i.e.,      */
            /* Transfer                                                  */
            /*                                                           */
            /*  p->displ=[ px1, px2;   to temp_m1 = pe = [ px1;          */
            /*             py1, py2;                       py1;          */
            /*             pz1, pz2 ]                      pz1;          */
            /*                                             px2;          */
            /*                                             py2;          */
            /*                                             pz2 ]         */
            /* --------------------------------------------------------- */

            temp_m1  = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 6, 1 );

            for( ii=0 ; ii < p->dof_per_node ; ++ii )
            for( jj=0 ; jj < p->nodes_per_elmt ; ++jj )
                 temp_m1->uMatrix.daa[ii+jj*p->dof_per_node][0] =
                    p->displ->uMatrix.daa[ii][jj];

            L = Element_Transformation_2d( p->coord, elmt_length );
            q = MatrixMult( L, temp_m1 );
            Q = MatrixMult( K, q );
            MatrixFree( temp_m1 );

            }

            if( p->elmt_state == 1 ) {
               L = Element_Transformation_2d( p->coord, elmt_length );
               Q = p->Q_saved;
            }