elmt_psps.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *          elmt_psps.c : Plane Stress Plane Strain 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: Mark Austin, Xiaoguang Chen, and Wane-Jang Lin      December 1995
 *  ============================================================================= 
 */

#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"

/*
#define DEBUG
*/
#define Streq(s1, s2) (strcmp(s1, s2) == 0)


/* ============================================================== */
/*   Element01                                                    */
/*   2D   Actions on Plane Stress/Plane Strain Element            */
/*        For Properties, see note below.                         */
/* ============================================================== */

save_actions_elmt01(ep,p)
ELEMENT *ep;
ARRAY    *p;
{
   /* ADD DETAILS LATER */
}


/* ================================================ */
/*   Plane-stress and plane strain Element          */
/*   elements: PLANE_STRAIN; PLANE_STRESS           */
/*   Plane Linear Elastic Element                   */
/* ================================================ */
/*    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;                                  */
/* ============================================================== */


ARRAY *elmt_psps(p,isw)
ARRAY *p;
int isw;
{
static QUANTITY fy, G, E, ET, density, lambda;
double nu;
double temperature, *alpha_thermal;

char *elmt_type;
int no_integ_pts, no_stress_pts;
int l, k, i, j,j1, k1, lint, kk;

double sg[17],tg[17],wg[17],sig[7],
       jacobain,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     *d1, *d2, *d3;

#ifdef DEBUG
       printf("*** enter elmt_psps()\n");
#endif

    UNITS_SWITCH = CheckUnits();
    dof_per_elmt  = p->nodes_per_elmt*p->dof_per_node;  

   /* Allocation for matrices */

    alpha_thermal = dVectorAlloc(2);
    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);
/*
   Mater_matrix  = MatrixAllocIndirectDouble(3, 3);
   temp_change   = MatrixAllocIndirectDouble(3, 1);
*/

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

            alpha_thermal[0] = p->work_material[7].value; /* thermal expansion coeff. in x-dirct */
            alpha_thermal[1] = p->work_material[8].value; /* thermal expansion coeff. in y-direc */

            if(p->nodes_per_elmt == 4) no_integ_pts  =  2;/* 4-node element */
            if(p->nodes_per_elmt == 8) no_integ_pts  =  3;/* 8-node element */

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

            elmt_type  = p->elmt_type;          

            if(Streq("PLANE_STRAIN", elmt_type)) {
               E.value  = p->work_material[0].value = E.value /(1+nu)/(1-nu);
               nu = p->work_material[4].value = nu * E.value;
               G.value  = p->work_material[1].value = E.value/2/(1+nu);
               for (i = 1; i <= 2; i++)
               alpha_thermal[i-1]  =  alpha_thermal[i-1]*(1.0 + nu);
            }

   switch(isw) {
       case PROPTY:  /* input material properties */

            lint = 0;
            break;
       case CHERROR:
            break;
       case STIFF:

#ifdef DEBUG
       printf("*** in elmt_psps() : start case STIFF \n");
#endif

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

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

            /* start gaussian integration  */

            for(l = 1; l<= lint; l++) { 
               shape(sg[l-1],tg[l-1],p->coord,shp,&jacobain,p->no_dimen,p->nodes_per_elmt, p->node_connect,FALSE);
            /* output:                                    */
            /*     shp[0][i-1] = dN_i/dx                  */
            /*     shp[1][i-1] = dN_i/dy                  */
            /*     compute [B] matrix at each Gaussian    */
            /*     integration point                      */

            /*****************************************************/
            /*  Derivative matrix                                */
            /*     (B_i)3x2; 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]          */
            /*       n = no of node                              */ 
            /*****************************************************/
            
            /* material matrix */
            Mater_matrix = MATER_MAT_PLANE(E.value, nu);

            /* 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];
            }
           
           /* muti. Jacobain determinant with weight coefficents and mater matrix */
          
            jacobain = jacobain* wg[l-1];

            for( i = 1; i <=3 ; i++ ) 
              for (j = 1; j <= 3; j++)  
                 Mater_matrix[i-1][j-1]  *= jacobain;

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

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

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

           /* [b] 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];
           }

#ifdef DEBUG
       dMatrixPrint("element stiffness matrix", p->stiff->uMatrix.daa, dof_per_elmt, dof_per_elmt);

      /* check element stiffness : sum of row of K = 0 */

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

       for (i = 1; i <= dof_per_elmt; i++) 
          for (j = 1; j <= dof_per_elmt; j++) 
            sum_row[i-1] +=  p->stiff->uMatrix.daa[i-1][j-1]; 

       for (i = 1; i <= dof_per_elmt; i++) 
         printf(" Stiffness K_e : sum of row[%d] = %lf \n", i, sum_row[i-1]); 

#endif

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

       /* Initiation of Stiffness Units Buffer                      */
     
       switch(UNITS_SWITCH) {
         case ON:
           if(UNITS_TYPE == 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]) );

           /* node i  i > 1*/
           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 OFF:
         break;
         default:
         break;
       }
     
#ifdef DEBUG
       printf("*** in elmt_psps() : leaving case STIFF\n");
#endif
       
      break;
      case EQUIV_NODAL_LOAD:
      /* calculate the equivalent nodal load */
      /* due to distributed loading          */
           
#ifdef DEBUG
       printf("*** in elmt_psps() : start case EQUIV_NODAL_LOAD: \n");
#endif
      /* initilize 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;
     }
      
     if(no_integ_pts*no_integ_pts != lint)
        pgauss(no_integ_pts,&lint,sg,tg,wg); 

     /* start gaussian integration */

     for(l = 1; l <= lint; l++) { 
         shape(sg[l-1],tg[l-1],p->coord,shp,&jacobain,p->no_dimen,p->nodes_per_elmt, p->node_connect,0);

            /* output: shp[0][i-1] = dN_i/dx          */
            /*         shp[1][i-1] = dN_i/dy          */
            /* compute [B] matrix at each Gaussian    */
            /* integration point                      */

            /*****************************************************/
            /*  Derivative matrix (B_i)3x2;                      */
            /*     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]          */
            /*      n  = no of node                              */ 
            /*****************************************************/
            
            /* material matrix */
            Mater_matrix = MATER_MAT_PLANE(E.value, nu);      /* need to be changed */
    
            /* 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];
            }
            
           /* muti. Jacobain determinant with weight coefficents and mater matrix */
            jacobain = jacobain* wg[l-1];
           
           /* [a] CALCULATE EQUIVALENT NODAL FORCE DUE TO INITIAL STRAIN     */

            if(p->nodal_init_strain != NULL) {

                for( i = 1; i <=3 ; i++ ) 
                  for (j = 1; j <= 3; j++)  
                    Mater_matrix[i-1][j-1]  *= jacobain;

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

            /* calculate strain[] at gaussian points : strain = sum Ni*nodal_strain */

                for (i = 1; i <= 3; i++)
                  strain[i-1][0]  =  0.0;
                
                for (i = 1; i <= 3; i++) {
                  for( j = 1; j <= p->nodes_per_elmt; j++) {
                    strain[i-1][0]  += shp[2][j-1] * p->nodal_init_strain[i-1][j-1];
                  }
                }
           /* [Mater]_3x3 * [strain]_3x1 and save in [stress]_3x1    */

                stress = dMatrixMultRep(stress, Mater_matrix, 3, 3, strain, 3, 1);
           
           /* mutiply [B]^T * [Mater]* [strain] */

              if(nodal_load == NULL ) { 
                  nodal_load = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
              } else {
                  load = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                  for (i = 1; i<= dof_per_elmt; i++)  {
                       nodal_load[i-1][0] += load[i-1][0];
                  }
              }
            }

           /* [b] CALCULATE EQUIVALENT NODAL FORCE DUE TO INITIAL STRESS     */

            if(p->nodal_init_stress != NULL) {

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

            /* calculate stress[] at gaussian points : stress = sum Ni*nodal_stress */

                for (i = 1; i <= 3; i++)
                  stress[i-1][0]  =  0.0;
                
                for (i = 1; i <= 3; i++) {
                  for( j = 1; j <= p->nodes_per_elmt; j++) {
                    stress[i-1][0]  += shp[2][j-1] * p->nodal_init_stress[i-1][j-1].value;
                  }
                  stress[i-1][0]  *=  jacobain;
                }

           /* mutiply [B]^T * [stress] */
               
                if(nodal_load == NULL) {
                  nodal_load  = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                } else {
                  load        = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                  for (i = 1; i<= dof_per_elmt; i++) 
                    nodal_load[i-1][0] -= load[i-1][0];
                  }
            }
           /* [c] CALCULATE EQUIVALENT NODAL FORCE DUE TO BODY FORCE         */

            if(p->nodal_body_force != NULL) {