elmt_shell_4n_q.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  elmt_shell_4n_q.c.
 *                                                                     
 *  This is a 4 node quadralateral shell element, incorportating membrane with
 *  normal drilling d.o.f. modified to remove effects of constant strains and
 *  including deep shell curvature corrections to the discrete kirchoff
 *  quadralateral plate 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 purpose, 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.
 *                                                                    
 *  ============================================================================= 
 *
 *   Element SHELL_FOUR_NODES
 *        Shell Element:
 *        material properties array
 *        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; 
 *
 *    p->no_dimen          = ndm;    ( = 3 ) (x,y,z cartesian coordinates at nodes)
 *    p->elmt_type         = &iel;   ( = 16 )
 *    p->nodes_per_elmt    = nen;    ( = 4 nodes counterclockwise around element)
 *    p->nodes_per_elmt    = nel;    ( = 4 )    = nen ???
 *    p->dof_per_node      = ndf;    ( = 6 )
 *    p->size_of_stiff     = nst;
 *    p->type              = d[16];   ???
 *    p->coord[ndm][nen]   = xl(ndm,*);
 *    p->node_connect[nen] = ix(nen);
 * 
 *    p->LC_ptr->ialph     = ialph;     Rotation constant
 *    p->LC_ptr->pen       = pen;       Penalty constant
 *    p->LC_ptr->load[0-5] = d[6-11];   Loading parameters
 *
 *    Young's Modulus = E.value
 *    Poisson Ratio   = nu
 *    Density         = density.value
 *    Shell thickness = h
 *          : 1-gravity       = d[6]
 *          : 2-gravity       = d[7]
 *          : pressure        = d[8]
 *          : x-gravity       = d[9]
 *          : y-gravity       = d[10]
 *          : z-gravity       = d[11]
 *          : Load Constant   = d[16]
 *          : Rotation Const. = (0=none, 1=hughes)\n", ialph
 *          : Penalty Const.  = pen
 *
 *  Written by: Lanheng Jin                                         December 1995
 *  Last modified Sept. 26, 1996               (Changes for stress output format)
 *  ============================================================================= 
 */

#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 */

static double b[4], c[4], aa[4], bb[4], cc[4], dd[4], ee[4];
static double bm[3][6][9];
static int    ii1, ii2;
static double shp[3][8][9], shp1[3][4][9], shp2[3][4][9];
static int    nlla, nllb, iel, ia[2][26], ir[2][26];
static int    false=0, true=1;
static int    mct = 0;

/*
 *  ======================================================================
 *  elmt_shell_4nodes_q() : 
 *  
 *  Input : 
 *  Output : 
 *  ======================================================================
 */

ARRAY *elmt_shell_4n_q(p, isw)
ARRAY *p;
int    isw;
{
static double                   nu, h; 
static QUANTITY                 G, E, density; 
static double *sg     = NULL;
static double *tg     = NULL;
static double *wg     = NULL;
static double **yl    = NULL;
static double **tr    = NULL;
static double *eps    = NULL;
static double *sigi   = NULL;
static double *sigb   = NULL;
static double **gshp1 = NULL;
static double **gshp2 = NULL;
static double **btd   = NULL;
static double **vl    = NULL;
static double *dvl    = NULL;
static double **s     = NULL;
static double *pp     = NULL;
static double *gg     = NULL;
static double *bl     = NULL;
static double *d      = NULL;

double        pen;
int           ialph;
int           l, lint;
int           qdflg;
int           ndm, ndf, nst, nen, nel;
int           i1, i2, j1, j2, ii, jj, kk;
int           i, j, k;
int           length1, length2, length;
double        a11, a12, a13, a21, a22, a23, a31, a32, a33, xx, yy, zz;
double        dv, dv1, dv2, dv3, dv4, dv5, dv6, ggi, ggv, tc, xsj, xsjt;
double        xyn, x1n, x2n, y1n, y2n, xn, yn, sn;
double        shp1i, shp2i, shp3i, shp11, shp12, shp13, shp21, shp22, shp23;
double        **ul;
DIMENSIONS    *dimen1, *dimen2, *dimen3;
double        coord, stress, moment, curvtr;

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

   /* [a] : Initialize element parameters */

   ndm = p->no_dimen;        /*    = 3      */ 
   nel = p->nodes_per_elmt;  /*    = 4      */
   nen = p->nodes_per_elmt;  /*    = 4      */
   ndf = p->dof_per_node;    /*    = 6      */
   nst = p->size_of_stiff;   /*    = nst    */

   if( isw != PROPTY && isw != CHERROR ) {

       l = d[14];
       if( l*l != lint )
          pgauss( l, &lint, sg, tg, wg );
       
       /* Compute the transformation and midsurface coordinates */

       if(yl == NULL)
          yl = MatrixAllocIndirectDouble(3,4);
       if(tr == NULL)
          tr = MatrixAllocIndirectDouble(3,3);

       tran06(p, yl, tr);

       /* Test for triangular element */
 
       if( p->node_connect[0] == p->node_connect[1] ||
           p->node_connect[1] == p->node_connect[2] ||
           p->node_connect[2] == p->node_connect[3] ||
           p->node_connect[3] == p->node_connect[0] ) {

          qdflg = false;
          for( i=1; i<=3; i++ ) {
             for( j=1; j<=6; j++ ) {
                for( k=1; k<=9; k++ ) {
                   bm[i-1][j-1][k-1] = 0.0;
                }
             }
          }

          shp11 = 0.0;
          shp12 = 0.0;
          shp13 = 0.0;
          shp21 = 0.0;
          shp22 = 0.0;
          shp23 = 0.0;
          a13   = 0.0;
          a23   = 0.0;
          a31   = 0.0;
          a32   = 0.0;
          a33   = 0.0;
          x1n   = 0.0;
          y1n   = 0.0;
          x2n   = 0.0;
          y2n   = 0.0;
          xyn   = 0.0;

          jtri06( yl, &xsjt );
       }
       else {
          qdflg = true;
          jacq06( yl );
       }

       if( qdflg && yl[2][0] == 0.0 ) {
          ii1 = 3;
          ii2 = 5;
       }
       else {
          ii1 = 1;
          ii2 = 6;
       }

       /* Construct integrals of the drilling shape functions */

       for( i=1; i<=3; i++ ) {
          for( j=1; j<=4; j++ ) {
             gshp1[i-1][j-1] = 0.0;
             gshp2[i-1][j-1] = 0.0;
          }
       }
       
       for( i=1; i<=24; i++ ) 
          gg[i-1] = 0.0;

       dv = 0.0;
     
       for( k=1; k<=lint; k++ ) {

          rshp06( k-1, sg[k-1], tg[k-1], yl, &xsj, ndm);

          dvl[k-1] = xsj * wg[k-1];
          dv      += dvl[k-1];
          for( i=1; i<=3; i++ ) {
             for( j=1; j<=4; j++ ) {
                gshp1[i-1][j-1]  += shp1[i-1][j-1][k-1]*dvl[k-1];
                gshp2[i-1][j-1]  += shp2[i-1][j-1][k-1]*dvl[k-1];
             }
          }
       }

       ggv = pen * d[2]/dv;
   }

   /*
    *  =======================================
    *  Beginning of the main switch 
    *  =======================================
    */

    switch(isw) {
        case PROPTY: 

             /* Allocate memory for working arrays */

             if(sg == NULL)
                sg = dVectorAlloc( 9 );
             if(tg == NULL)
                tg = dVectorAlloc( 9 );
             if(wg == NULL)
                wg = dVectorAlloc( 9 );
             if(gg == NULL)
                gg = dVectorAlloc( 24 );
             if(dvl == NULL)
                dvl = dVectorAlloc( 9 );
             if(pp == NULL)
                pp = dVectorAlloc(nst);
             if(gshp1 == NULL)
                gshp1 = MatrixAllocIndirectDouble( 3, 4 );
             if(gshp2 == NULL)
                gshp2 = MatrixAllocIndirectDouble( 3, 4 );
             if(btd == NULL) 
                btd = MatrixAllocIndirectDouble(3,6);
             if(s == NULL)
                s = MatrixAllocIndirectDouble(nst, nst);

             /* Material properties: elastic */

             E.value       =  p->work_material[0].value;
             nu            =  p->work_material[4].value;
             density.value =  p->work_material[5].value;

             /* Check value of Poisson's ratio */

             if( nu == 0.0 || nu > 0.5 ) {
                 printf("WARNING >> ... In elmt_shell_4nodes_q() : nu value = %9.4f\n", nu);
                 printf("WARNING >> ... Reset nu to 0.3 !\n");
                 nu = 0.3;       
             }

             /* Calculate value of G */

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

             if( CheckUnits() == ON ) {
                 E.dimen       =  p->work_material[0].dimen;
                 density.dimen =  p->work_material[5].dimen;
                 G.dimen = E.dimen;
             }

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

             /* Get thickness of the shell  */

             h = p->work_section[11].value;  

             if(d == NULL)
                d = dVectorAlloc( 20 );

             ialph = p->LC_ptr->ialph;
             pen   = p->LC_ptr->pen;
             d[6]  = p->LC_ptr->load[0];
             d[7]  = p->LC_ptr->load[1];
             d[8]  = p->LC_ptr->load[2];
             d[9]  = p->LC_ptr->load[3];
             d[10] = p->LC_ptr->load[4];
             d[11] = p->LC_ptr->load[5];
             d[0]  = E.value/(1.0-nu*nu)*h;
             d[1]  = nu*d[0];
             d[2]  = (d[0] - d[1])/2.0;

             xx    = h*h/12.0;

             d[3]  = d[0]*xx;
             d[4]  = d[1]*xx;
             d[5]  = d[2]*xx;
             d[12] = density.value*h;
             d[13] = d[12]*xx;

             /* Quadrature order */

             d[14] = 2.0;
             d[15] = 1.0;
             lint  = 0;

             /* Construct rotation parameters : u-x = 1, u-y = 2 */

             iel = 16;       /* set up the element type 16   */

             ia[0][iel-1] = 1;
             ia[1][iel-1] = 2;

             /* Construct rotation parameters : theta-x = 4, theta-y = 5 */

             ir[0][iel-1] = 4;
             ir[1][iel-1] = 5;

             break;

        case CHERROR:
             break;
        case STIFF:         /* Compute the element tangent stiffness array  */
        case LOAD_MATRIX:   /* Form load matrix */

             /* Construct the modified drilling shape functions */

             for( i=1; i<=3; i++ ) {
             for( j=1; j<=4; j++ ) {
                  dv1 = gshp1[i-1][j-1]/dv;
                  dv2 = gshp2[i-1][j-1]/dv;
                  for( k=1; k<=lint; k++ ) {
                     shp1[i-1][j-1][k-1] -= dv1;
                     shp2[i-1][j-1][k-1] -= dv2;
                  }
             }
             }

             /* Transform the element load vector (dm is proportional to load factor) */

             if(bl == NULL)
                bl = dVectorAlloc(3);

             for( i=1; i<=3; i++ ) {
                bl[i-1] = d[i+5];                   /*  set dm=1 here  */
                for( j=1; j<=3; j++ )
                   bl[i-1]  += tr[i-1][j-1]*d[i+8]; /*  set dm=1 here  */
             }

             tc = d[17-1];

             /* Compute membrane/load parts */
 
             for( k=1; k<=lint; k++ ) {

                /* Membrane and bending stiffness parts */

                if( qdflg ) 
                    dktq06( k-1 );
                else 
                    dktb06( sg[k-1], tg[k-1], xsjt );
         
                dv  = dvl[k-1];
                dv1 = d[1-1] * dv;
                dv2 = d[2-1] * dv;
                dv3 = d[3-1] * dv;
                dv4 = d[4-1] * dv;