elmt_fiber3d.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  elmt_fiber3d.c : Linear/Nonlinear 3D 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.
 *                                                                    
 *  Written by: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin         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
*/

#ifdef __STDC__
ARRAY *elmt_fiber_3d(ARRAY *p, int isw)
#else
ARRAY *elmt_fiber_3d(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      cx, cy, cz, cxz, xx, yy, zz, temp1, temp2, temp3;
double      **rot, mbar, rT;
MATRIX      *temp_m1, *temp_m2;
DIMENSIONS  *dp_length, *dp_force, *dp_stress, *dimen;

double      G, J, bf, depth;
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();

    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;
          cx = p->coord[0][1].value - p->coord[0][0].value;
          cy = p->coord[1][1].value - p->coord[1][0].value;
          cz = p->coord[2][1].value - p->coord[2][0].value;
          elmt_length.value = sqrt( cx*cx + cy*cy + cz*cz );
          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 );

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

          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_3d( bx, xi, elmt_length );
             Section_Tangent_Stiffness_3d( 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 );
          }
          K = MatrixInverse( F );

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

          /* Rigid body rotation transformation */
          /* Ke = Trans(R)*K*R , element stiffness in local coordinate system */
          R = Rigid_Body_Rotation_3d( elmt_length );
          Rtrans = MatrixTranspose( R );
          temp_m1 = MatrixMult( Rtrans, K );
          Ke = MatrixMult( temp_m1, R );
          MatrixFree( temp_m1 );

          /* Calculate torsional stiffness */
          /* Assuming plane remain plane, so twisting stiffness is independent to others */
          if ( (G=p->work_material[1].value) <= 0.0 )
             G = p->work_material[0].value/2.0/(1+p->work_material[4].value);
          J = p->work_section[12].value;
          bf = p->work_section[7].value;
          depth = p->work_section[9].value;

          /* If J value not input, J calculated based on rectangular */
          /* section of size (bf x depth)                            */

          if(J == 0.0 ) {
             if(bf == 0.0 || depth == 0.0){
                printf("WARNING >> Must give 'J' or ('width' & 'depth') to calculate stiffness");
                exit(1);
             }
             /* Check bf < depth & cal J */
             if(bf < depth ) 
                J = (1.0-0.63*bf/depth)*bf*bf*bf*depth/3.0;
             else
                J = (1.0-0.63*depth/bf)*depth*depth*depth*bf/3.0;
          }

          Ke->uMatrix.daa[3][3] = G*J/elmt_length.value;
          Ke->uMatrix.daa[3][9] = -Ke->uMatrix.daa[3][3];
          Ke->uMatrix.daa[9][3] =  Ke->uMatrix.daa[3][9];
          Ke->uMatrix.daa[9][9] =  Ke->uMatrix.daa[3][3];

          /* Transform local coordinate system to global coordinate system */

          cx = (p->coord[0][1].value - p->coord[0][0].value)/elmt_length.value;
          cy = (p->coord[1][1].value - p->coord[1][0].value)/elmt_length.value;
          cz = (p->coord[2][1].value - p->coord[2][0].value)/elmt_length.value;
          cxz = sqrt( cx*cx + cz*cz );
          if( cx != 1.0 ) {
             if( cxz == 0.0 ) {
                /* Ke' = Ke*T */
                for( ii=0 ; ii < 12 ; ++ii ) {
                   for( jj=0 ; jj < 12 ; jj=jj+3 ) {
                      temp1 = Ke->uMatrix.daa[ii][jj];
                      temp2 = Ke->uMatrix.daa[ii][jj+1];
                      Ke->uMatrix.daa[ii][jj]   = -temp2*cy;
                      Ke->uMatrix.daa[ii][jj+1] =  temp1*cy;
                   }
                }
                /* Ke(global) = Trans(T)*Ke*T = Trans(T)*Ke' */
                for( ii=0 ; ii < 12 ; ii=ii+3 ) {
                   for( jj=0 ; jj < 12 ; ++jj ) {
                      temp1 = Ke->uMatrix.daa[ii][jj];
                      temp2 = Ke->uMatrix.daa[ii+1][jj];
                      Ke->uMatrix.daa[ii][jj]   = -temp2*cy;
                      Ke->uMatrix.daa[ii+1][jj] =  temp1*cy;
                   }
                }
             }
             else {
                /* Ke' = Ke*T */
                for( ii=0 ; ii < 12 ; ++ii ) {
                   for( jj=0 ; jj < 12 ; jj=jj+3 ) {
                      temp1 = Ke->uMatrix.daa[ii][jj];
                      temp2 = Ke->uMatrix.daa[ii][jj+1];
                      temp3 = Ke->uMatrix.daa[ii][jj+2];
                      Ke->uMatrix.daa[ii][jj]   = temp1*cx - temp2*cx*cy/cxz - temp3*cz/cxz;
                      Ke->uMatrix.daa[ii][jj+1] = temp1*cy + temp2*cxz;
                      Ke->uMatrix.daa[ii][jj+2] = temp1*cz - temp2*cy*cz/cxz + temp3*cx/cxz;
                   }
                }
                /* Ke(global) = Trans(T)*Ke*T = Trans(T)*Ke' */
                for( ii=0 ; ii < 12 ; ii=ii+3 ) {
                   for( jj=0 ; jj < 12 ; ++jj ) {
                      temp1 = Ke->uMatrix.daa[ii][jj];
                      temp2 = Ke->uMatrix.daa[ii+1][jj];
                      temp3 = Ke->uMatrix.daa[ii+2][jj];
                      Ke->uMatrix.daa[ii][jj]   = temp1*cx - temp2*cx*cy/cxz - temp3*cz/cxz;
                      Ke->uMatrix.daa[ii+1][jj] = temp1*cy + temp2*cxz;
                      Ke->uMatrix.daa[ii+2][jj] = temp1*cz - temp2*cy*cz/cxz + temp3*cx/cxz;
                   }
                }
             }
          }

          /* 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( K );
          MatrixFree( Rtrans );
          MatrixFree( R );
          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]) );
             ZeroUnits( &(p->stiff->spRowUnits[2]) );
             UnitsCopy( &(p->stiff->spRowUnits[3]), dp_length );
             UnitsCopy( &(p->stiff->spRowUnits[4]), dp_length );
             UnitsCopy( &(p->stiff->spRowUnits[5]), dp_length );
             UnitsCopy( &(p->stiff->spRowUnits[6]),  &(p->stiff->spRowUnits[0]) );
             UnitsCopy( &(p->stiff->spRowUnits[7]),  &(p->stiff->spRowUnits[1]) );
             UnitsCopy( &(p->stiff->spRowUnits[8]),  &(p->stiff->spRowUnits[2]) );
             UnitsCopy( &(p->stiff->spRowUnits[9]),  &(p->stiff->spRowUnits[3]) );
             UnitsCopy( &(p->stiff->spRowUnits[10]), &(p->stiff->spRowUnits[4]) );
             UnitsCopy( &(p->stiff->spRowUnits[11]), &(p->stiff->spRowUnits[5]) );

             UnitsDivRep( &(p->stiff->spColUnits[0]), dp_force, dp_length, NO );
             UnitsCopy(   &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]) );
             UnitsCopy(   &(p->stiff->spColUnits[2]), &(p->stiff->spColUnits[0]) );
             UnitsCopy(   &(p->stiff->spColUnits[3]), dp_force );
             UnitsCopy(   &(p->stiff->spColUnits[4]), dp_force );
             UnitsCopy(   &(p->stiff->spColUnits[5]), dp_force );
             UnitsCopy( &(p->stiff->spColUnits[6]),  &(p->stiff->spColUnits[0]) );
             UnitsCopy( &(p->stiff->spColUnits[7]),  &(p->stiff->spColUnits[1]) );
             UnitsCopy( &(p->stiff->spColUnits[8]),  &(p->stiff->spColUnits[2]) );
             UnitsCopy( &(p->stiff->spColUnits[9]),  &(p->stiff->spColUnits[3]) );
             UnitsCopy( &(p->stiff->spColUnits[10]), &(p->stiff->spColUnits[4]) );
             UnitsCopy( &(p->stiff->spColUnits[11]), &(p->stiff->spColUnits[5]) );

             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:
       case STRESS:

          if( UNITS_SWITCH == ON )
             SetUnitsOff();

          if( p->elmt_state == 0 ) {  /* Never pass the ElmtStateDet() */

             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 );

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

             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_3d( bx, xi, elmt_length );
                Section_Tangent_Stiffness_3d( 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 );
             }
             K = MatrixInverse( F );

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