elmt_psps.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  elmt_psps.c : Plane Stress Plane Strain 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 "units.h"
#include "matrix.h"
#include "vector.h"
#include "fe_database.h"
#include "symbol.h"
#include "fe_functions.h"
#include "elmt.h"

/* Function Declarations */

double **MaterialMatrixPlane( double, double, double );

/*
 *  ============================================================================ 
 *  Plane-stress/plane strain finite element         
 * 
 *  The material and section properties are: 
 * 
 *      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_material[7] = alpha_thermal[0];
 *      p->work_material[8] = alpha_thermal[1];
 *
 *      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;                                 
 * 
 *  Compute (3x2) derivative matrix [B] .... 
 * 
 *       B_i[0][0] = dNi/dx, B_i[0][1] = 0  
 *       B_i[1][0] = 0,      B_i[1][1] = dNi/dy     
 *       B_i[2][0] = dNi/dy, B_i[2][2] = dNi/dx   
 *
 *       [B] = [B_1, B_2, B_3, B_4, ..., B_n] where n = no of node                             
 *
 *       The output is : shp[0][i-1] = dN_i/dx                
 *                       shp[1][i-1] = dN_i/dy                
 *       compute [B] matrix at each Gaussian integration point.
 *
 *  Note : Material properties such and "E" and "nu" must be stored as static
 *         variables.
 *  ============================================================================ 
 */

ARRAY *elmt_psps(p,isw)
ARRAY *p;
int isw;
{
static QUANTITY fy, G, E, ET, density, lambda;
static double nu;
static double temperature;
static double *alpha_thermal;
static double dFlag;
static int no_integ_pts;
static int no_stress_pts;
int ii, k, i, j,j1, k1, lint, kk;
double sg[17],tg[17],wg[17],sig[7],
       jacobian,w11,w12,w21,w22, xx, yy, dv, wd[3];
double **B_matrix;
double **B_Transpose;
double **stiff;
double **Mater_matrix;
double **body_force;
double **strain;
double **stress;
double **temp_change;
double **load;
double **displ;
double **nodal_load;
double **shp;
double **temp_matrix;

int             dof_per_elmt;
int          length1, length;
int             UNITS_SWITCH;
double         *sum_row, sum;
DIMENSIONS        *dp_stress;
DIMENSIONS        *dp_length;
DIMENSIONS     *d1, *d2, *d3;

   dof_per_elmt  = p->nodes_per_elmt*p->dof_per_node;  

    /* [a] : Allocate memory for matrices */

   strain        = MatrixAllocIndirectDouble(3, 1);
   stress        = MatrixAllocIndirectDouble(3, 1);
   displ         = MatrixAllocIndirectDouble(dof_per_elmt, 1);
   stiff         = MatrixAllocIndirectDouble(dof_per_elmt, dof_per_elmt);
   load          = MatrixAllocIndirectDouble(dof_per_elmt, 1);
   nodal_load    = MatrixAllocIndirectDouble(dof_per_elmt, 1);
   B_matrix      = MatrixAllocIndirectDouble(3, dof_per_elmt);
   shp           = MatrixAllocIndirectDouble(3, p->nodes_per_elmt);
   sum_row       = dVectorAlloc(dof_per_elmt);

   B_Transpose   = MatrixAllocIndirectDouble(dof_per_elmt, 3);
   temp_matrix   = MatrixAllocIndirectDouble(dof_per_elmt, 3);
   body_force    = MatrixAllocIndirectDouble(3, 1);

   /* [b] : deal with Individual Tasks */

   UNITS_SWITCH = CheckUnits();
   switch(isw) {
       case PROPTY: 
            alpha_thermal = dVectorAlloc(2);

            /* Input material properties */

            lint = 0;
            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;
            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;
            }

            /* Initialize thermal expansion coefficients in x- and y- directions */

            alpha_thermal[0] = p->work_material[7].value; 
            alpha_thermal[1] = p->work_material[8].value; 

            /* 4-node and 8-node elements */

            if(p->nodes_per_elmt == 4)
               no_integ_pts  =  2;
            if(p->nodes_per_elmt == 8)
               no_integ_pts  =  3;

            no_stress_pts = no_integ_pts; /* stress points, subject to change later */

            /* Set flag for Plane Stress/Plane Strain material matrix computations */

            dFlag = 1;
            if(strcmp("PLANE_STRAIN", p->elmt_type) == 0) {
               dFlag = 2;
            }

            break;
       case CHERROR:
            break;
       case STIFF: 

            /* Zero out the Stiffness Matrix */

            for(i = 1; i <= dof_per_elmt; i++) 
                for(j = 1; j<= dof_per_elmt; j++)
                    stiff[i-1][j-1] = 0.0;

            /* Get Gauss Integration points */

            if(no_integ_pts*no_integ_pts != lint)
               pgauss(no_integ_pts,&lint,sg,tg,wg);

            /* Start Gaussian Integration  */

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

                /* Compute shape functions */

                shape( sg[ii-1],tg[ii-1], p->coord,shp,&jacobian, p->no_dimen,
                       p->nodes_per_elmt, p->node_connect,FALSE);

                /*
                 *  =================================================== 
                 *  Compute material matrix
                 *  
                 *  Plane Stress : 3rd argument ii = 1
                 *  Plane Strain : 3rd argument ii = 2
                 *  =================================================== 
                 */

                Mater_matrix = MaterialMatrixPlane( E.value, nu, dFlag );

                /* Form [B] matrix */
            
                for(j = 1; j <= p->nodes_per_elmt; j++) { 
                    k = 2*(j-1);
                    B_matrix[0][k]   = shp[0][j-1];
                    B_matrix[1][k+1] = shp[1][j-1];
                    B_matrix[2][k]   = shp[1][j-1];
                    B_matrix[2][k+1] = shp[0][j-1];
                }
           
                /* Multiply Jacobian determinant with weight coefficents  */
                /* and material matrix                                    */
          
                jacobian = jacobian*wg[ii-1];
                for( i = 1; i <=3 ; i++ ) 
                    for (j = 1; j <= 3; j++)  
                         Mater_matrix[i-1][j-1] *= jacobian;

                /* Transpose [B] matrix */

                for(i = 1; i <= 3; i++)
                    for(j = 1; j <= dof_per_elmt; j++) 
                        B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];

                /* Compute [B]^T*[Mater] and save in [B]^T */

                temp_matrix = dMatrixMultRep(temp_matrix, B_Transpose,
                                             dof_per_elmt, 3, Mater_matrix, 3, 3);

                /* Calculate stiffness : Integral of [B]^T * [Mater] * [B] */

                stiff = dMatrixMultRep(stiff, temp_matrix, dof_per_elmt, 3,
                                       B_matrix, 3, dof_per_elmt);
 
                for (i = 1; i <= dof_per_elmt; i++) 
                for (j = 1; j <= dof_per_elmt; j++) 
                     p->stiff->uMatrix.daa[i-1][j-1]  += stiff[i-1][j-1];

                /* Free memory for material matrix with modified contents */

                MatrixFreeIndirectDouble(Mater_matrix, 3);
            }

            /* [c] : Assign units to stiffness matrix */
     
            if ( UNITS_SWITCH == ON ) {
               if( CheckUnitsType() == SI) {
                   d1 = DefaultUnits("Pa");
                   d2 = DefaultUnits("m");
               } else {
                   d1 = DefaultUnits("psi");
                   d2 = DefaultUnits("in");
               }

               /* Node 1 */

               UnitsMultRep( &(p->stiff->spColUnits[0]), d1, d2 );
               UnitsCopy( &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]) );

               /* Nodes 2 through 4 */

               for(i = 2; i <= p->nodes_per_elmt; i++) {
               for(j = 1; j <= p->dof_per_node; j++) {
                   k  = p->dof_per_node*(i-1) + j;
                   UnitsCopy( &(p->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) );
               }
               }

               free((char *) d1->units_name);
               free((char *) d1);
               free((char *) d2->units_name);
               free((char *) d2);

               for( i=1 ; i<=p->size_of_stiff ; i++ )
                    ZeroUnits( &(p->stiff->spRowUnits[i-1]) );

            }
            break;
       case EQUIV_NODAL_LOAD: /* calculate the equivalent nodal loads */
                              /* due to distributed loadings          */
           
            /* Initialize nodal_load */
      
            for(i = 1; i <= dof_per_elmt; i++) {
                p->equiv_nodal_load->uMatrix.daa[i-1][0] = 0.0;
                nodal_load[i-1][0] = 0.0;
                load[i-1][0]       = 0.0;
            }
 
            for(i = 1; i<= no_integ_pts*no_integ_pts; i++) {
                sg[i-1] = 0.0; tg[i-1] = 0.0; wg[i-1] = 0.0;
            }

            /* Compute Gaussian Integration Points */
      
            if(no_integ_pts*no_integ_pts != lint)
               pgauss(no_integ_pts,&lint,sg,tg,wg); 

            /* Main loop for Gauss Integration */

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

                /* Compute shape functions */

                shape(sg[ii-1],tg[ii-1],p->coord,shp,&jacobian,
                      p->no_dimen,p->nodes_per_elmt, p->node_connect,0);

                /* Form [B] matrix */
            
                for(j = 1; j <= p->nodes_per_elmt; j++) { 
                    k = 2*(j-1);
                    B_matrix[0][k]   = shp[0][j-1];
                    B_matrix[1][k+1] = shp[1][j-1];
                    B_matrix[2][k]   = shp[1][j-1];
                    B_matrix[2][k+1] = shp[0][j-1];
                }

                /* Get material matrix : see notes on dFlag above */

                Mater_matrix = MaterialMatrixPlane( E.value, nu, dFlag );