elmt_shell_4n_q.c

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

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *    elmt_shell_4n_q.c : 4 nodes shell quadratic 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.
 *                                                                    
 *   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 )
 *    p->elmt_type         = &iel;   ( = 16 )
 *    p->nodes_per_elmt    = nen;    ( = 4 )
 *    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
 *
 *  Written by: Lanheng Jin                                         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 
*/

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;

ARRAY *elmt_shell_4nodes_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;

FILE *fopen(), *fp;

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

    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 );
       
/* -- Calculation of transformation array and surface coordiantes -- */

       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 the 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;
    }

/**********************************************/
/**    Begining of the main switch           **/
/**********************************************/

    switch(isw) {

        case PROPTY: 

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

        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 possion_ratio value ------------------ */

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

/* ---- calculate  G value ------------------------- */

        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), check G for homogeneous material \n");
           printf("               Ignore this message for non-homogeneous materials \n");
        }

        /* Get the thickness of the shell  */

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

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

/*
        fp = fopen("data", "r");

        if(fp) fscanf(fp, "%d%lf%lf%lf%lf%lf%lf%lf%lf", 
                      &ialph, &pen, d+6, d+7, d+8, d+9, d+10, d+11, d+16);
        else {
           printf("\nError for opening a file\n");
           exit(1);
        }

        fclose (fp);
*/

	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;

/* ---- Print out some parameters ---------------------- */
/*
        printf("Output of material properties, section parameters, load and some constants\n");
        printf("    : Young's Modulus =   %l3.5e\n", E.value);
        printf("    : Poisson Ratio   = %13.5e\n", nu);
        printf("    : Density         = %13.5e\n", density.value);
        printf("    : Shell thickness = %13.5e\n", h);
        printf("    : 1-gravity       = %13.5e\n", d[6]);
        printf("    : 2-gravity       = %13.5e\n", d[7]);
        printf("    : pressure        = %13.5e\n", d[8]);
        printf("    : x-gravity       = %13.5e\n", d[9]);
        printf("    : y-gravity       = %13.5e\n", d[10]);
        printf("    : z-gravity       = %13.5e\n", d[11]);
        printf("    : Load Constant   = %13.5e   (0=lump, 1=cons.)\n", d[16]);
        printf("    : Rotation Const. = %3d             (0=none, 1=hughes)\n", ialph);
        printf("    : Penalty Const.  = %13.5e\n\n", pen);
*/
#ifdef DEBUG
    printf(" In elmt_shell_4nodes_q(): \n");
    printf("    : leaving case of PROPTY\n");
#endif
        break;

        case CHERROR:
        break;

        case STIFF:              /* form element stiffness */
                                 /* Compute the element tangent array  */
        case LOAD_MATRIX:        /* form load matrix */

#ifdef DEBUG
    printf("*** elmt_shell_4nodes_q() : In case STIFF\n");
#endif

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

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