elmt_shell_8n.c

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

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *      elmt_shell_8n.c : Eight Node Shell Finite Element
 *                                                                     
 *  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: Xiaoguang Chen                                      December 1995
 *  ============================================================================= 
 */

#include <stdio.h>
#include <string.h>
#include <math.h>

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

/*
#define DEBUG 
#define DEBUG1 
*/


/* ============================================================== */
/*   Element SHELL_EIGHT_NODES                                    */
/*        Implicit code                                           */
/*        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] = thickness;                            */
/* ============================================================== */

ARRAY *elmt_shell_8n(p, isw)
ARRAY    *p;
int     isw;
{
#ifdef DEBUG
       printf(" enter elmt_shell_8n() \n");
#endif

   p = elmt_shell_8nodes_implicit(p, isw);

#ifdef DEBUG
       printf(" Leaving elmt_shell_8n() \n");
#endif

   return(p);
}

ARRAY *elmt_shell_8nodes_implicit(p, isw)
ARRAY    *p;
int     isw;
{
static double                                   nu; 
static QUANTITY      fy, G, E, ET, density, lambda; 
static double         Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
static double        bf, tf, A, depth, h, EA, EIzz;

static DIMENSIONS *dp_length, *dp_force, *dp_moment;
static DIMENSIONS *dp_stress, *dp_degree, *dp_temperature;

double             EE, jacobian, k_bar = 0.833;     /* k_bar is shear factor */
double                   elmt_length, aspect_ratio;

double                           **co_coord = NULL;     
double                          **Direction_Matrix;
double                   **T_matrix, **T_Transpose;
double                             **e1_ptr = NULL;
double                             **e2_ptr = NULL;
double                             **e3_ptr = NULL;

static double      *integ_coord = NULL, *weight = NULL;

double                 **stress = NULL, **displ = NULL;
double   **nodal_load = NULL, **nodal_load_temp = NULL;
double                                 diff, sum, temp;
double     **Cep = NULL, **stiff = NULL, **mass = NULL; 
int             i, j, k, ii, jj, kk, n, n1, n2, k1, k2;
int          surface_pts, surf_integ_pts, no_integ_pts;
int                size, dof, length, length1, length2;
int                                       UNITS_SWITCH;


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

     
       UNITS_SWITCH = CheckUnits();
       co_coord         = MatrixAllocIndirectDouble(p->no_dimen, p->nodes_per_elmt);

       for( i = 1; i <= 3; i++ ) {
           for ( j = 1; j <= p->nodes_per_elmt; j++) {
               co_coord[i-1][j-1] = p->coord[i-1][j-1].value;
           }
       }
#ifdef DEBUG
        dMatrixPrint("co_coord", co_coord, 3,8);
#endif

      no_integ_pts = p->integ_ptr->thickness_pts;
      dof          = p->dof_per_node;
      size         = p->size_of_stiff;

/*****************************************************************/
#ifdef DEBUG
    printf(" ***** In elmt_shell_8nodes_implicit(): \n");
    printf("                                      : enter main switch \n");
#endif

    switch(isw) {
        case PROPTY: 

#ifdef DEBUG
    printf(" In elmt_shell_8nodes_implicit(): \n");
    printf("    : enter case of PROPTY\n");
#endif

             /* material properties:  elastic */

             E.value       =  p->work_material[0].value;
             fy.value      =  p->work_material[2].value;
             ET.value      =  p->work_material[3].value;
             nu            =  p->work_material[4].value;
             density.value =  p->work_material[5].value;
             lambda.value  =  nu*E.value/(1.0+nu)/(1.0-2.0*nu);

           /* (1)   check  possion_ratio value */

            if( nu == 0.0 || nu > 0.5 ) {
                printf("WARNING >> ... In shell_Belytschko element() -  nu value = %9.4f,reset to 0.3 !\n",
                       nu);
                nu = 0.3;    /* default poi_ratio value */
            }

            /* (2)   calculate  G value */
            
            G.value = p->work_material[1].value = E.value/(1.0 - 2.0*nu) ;

            if(E.value/((1.0 - 2.0*nu)) != p->work_material[1].value) {
                    printf(" elmt_shell(): WARNING: G is not equal to E/(1-2nu), check G for homogeneous material \n");
                    printf("             : ignore this message for non-homogeneous materials \n");
            }


            if( UNITS_SWITCH == ON ) {
               E.dimen       =  p->work_material[0].dimen;
               fy.dimen      =  p->work_material[2].dimen;
               ET.dimen      =  p->work_material[3].dimen;
               density.dimen =  p->work_material[5].dimen;
               lambda.dimen  =  E.dimen;
               G.dimen = E.dimen;
            }

             Ixx    = p->work_section[0].value;
             Iyy    = p->work_section[1].value;
             Izz    = p->work_section[2].value;
             Ixy    = p->work_section[3].value;
             Ixz    = p->work_section[4].value;
             Iyz    = p->work_section[5].value;
             bf     = p->work_section[7].value;
             tf     = p->work_section[8].value;
             depth  = p->work_section[9].value;
             A      = p->work_section[10].value;
             h      = p->work_section[11].value;    /* thickness of the shell */

             EA     = E.value*A;
             EIzz   = E.value*Izz;
         
#ifdef DEBUG
    printf(" In elmt_shell_8nodes_implicit(): \n");
    printf("    : leaving case of PROPTY\n");
#endif
             break;
             case STIFF: /* form element stiffness */

#ifdef DEBUG
       printf("*** In elmt_shell() : start case STIFF\n");
       printf("                    : Density         = %14.4e\n", density.value);
       printf("                    : shell_thickness = %8.2f\n", h);
       printf("                    : E               = %lf\n", E.value);
       printf("                    : nu              = %lf\n", nu);
       printf("                    : no_integ_pts in z_direction = %d\n", no_integ_pts);
#endif
              /* ============================================*/
              /* initialize the stiffness matrix             */
              /* ============================================*/
                 for(i = 1; i <= p->size_of_stiff; i++) {
                    for(j = 1; j <= p->size_of_stiff; j++) {
                        p->stiff->uMatrix.daa[i-1][j-1] = 0.0;
                    }
                 }

                 stiff        = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff);
                 integ_coord  = dVectorAlloc(no_integ_pts + 1);
                 weight       = dVectorAlloc(no_integ_pts + 1);

                 size         = p->size_of_stiff;
                 gauss(integ_coord,weight,no_integ_pts);

                 /* Integration loop over thickness */

                 for (ii = 1; ii <= no_integ_pts; ii++) {

                      Shell_Stiff_Plane_8node(stiff, p, co_coord, integ_coord[ii], ii, E.value, nu); 	

                      for (i = 1; i <= size; i++) {
                          for (j = 1; j <= size; j++) {
                               p->stiff->uMatrix.daa[i-1][j-1] += stiff[i-1][j-1]*weight[ii];     
                          }
                      }
                 }  /* end of gaussian integration through thickness */

                 free((char *) integ_coord);
                 free((char *) weight);

                MatrixFreeIndirectDouble(stiff, size);

                /**************************************************/
                /* Assign Units to Stiffness Matrix               */
                /**************************************************/

       /* Initiation of Stiffness Units Buffer                      */

       switch( UNITS_SWITCH ) {
         case ON:
           if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
              dp_stress = DefaultUnits("Pa");
              dp_length = DefaultUnits("m");
           }
           else {
              dp_stress = DefaultUnits("psi");
              dp_length = DefaultUnits("in");
           }

          /* node no 1 */
           UnitsMultRep( &(p->stiff->spColUnits[0]), dp_stress, dp_length );
           UnitsCopy( &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]) );
           UnitsCopy( &(p->stiff->spColUnits[2]), &(p->stiff->spColUnits[0]) );
           UnitsMultRep( &(p->stiff->spColUnits[3]), &(p->stiff->spColUnits[0]), dp_length );
           UnitsCopy( &(p->stiff->spColUnits[4]), &(p->stiff->spColUnits[3]) );

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

          /* node i  i > 1*/
        
           for ( i = 2; i <= p->nodes_per_elmt; i++) {
                kk = p->dof_per_node*(i-1) + 3; 
                for( j = 1; j <= p->dof_per_node; j++) {
                     k  = p->dof_per_node*(i-1) + j;
                     if( k <= kk) {
                        UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) );
                        UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[0]) );
                     }
                     if(k > kk) {
                        UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[3]) );
                        UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[3]) );
                     }
                 }
            }
            free((char *) dp_stress->units_name);
            free((char *) dp_stress);
            free((char *) dp_length->units_name);
            free((char *) dp_length);
          break;
          case OFF:
          break;
          default:
          break;
        }
             break;
	case PRESSLD:
	case EQUIV_NODAL_LOAD:

#ifdef DEBUG
       printf("*** In elmt_shell() : enter case EQUIV_NODAL_LOAD\n");
#endif
         /* Form external nodal load vector */
         /* due to distributed loading      */
            
       if(p->elmt_load_ptr != (ELEMENT_LOADS *) NULL) {
           p = sld108(p, PRESSLD); /* Equivalent Load in local_coordinate */
       }

/* ------------------------ UNITS -----------------------------------*/
      UNITS_SWITCH = CheckUnits();
      switch( UNITS_SWITCH ) {
        case ON:
         if(UNITS_TYPE == SI) {
            dp_length = DefaultUnits("m");
            dp_force  = DefaultUnits("N");
         }
         if(UNITS_TYPE == US) {
            dp_length = DefaultUnits("in");
            dp_force  = DefaultUnits("lbf");
         }
         dp_moment = UnitsMult( dp_force, dp_length );
             
         for(i= 1; i<= p->dof_per_node; i++) {
            for(j = 1; j <= p->nodes_per_elmt; j++) {
                k  = p->dof_per_node*(j-1)+i;
                k1 =  p->dof_per_node*j-3;
                if ( k <= k1) {      /* force units */
                   UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), dp_force );
                }
                else {  /* k > k1 moment units */
                   UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), dp_moment );
                }
             }
          }