fe_nonlinear.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  fe_nonlinear.c : Functions Needed for Nonlinear Finite Element Analysis
 *                                                                     
 *  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           May 1997
 *  Modified by: Mark Austin                                           March 2000
 *  ============================================================================= 
 */

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

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

extern ARRAY     *array;
extern EFRAME    *frame;

/* #define DEBUG */

/* ======================================================================= */
/* Setup pointers to working arrays : notes static flag limits             */
/* scope of arrays to fe_nonlinear.c (i.e., don't break this file apart!). */
/* ======================================================================= */

static RESPONSE                *RespondBuff;
static int                   *ElmtStateBuff;
static MATER_LOAD_CURVE          *LoadCurve;
static FIBER_RESPONSE   *FiberRespondBuffer;
static int                  **FiberLoadFlag;

/*
 * =======================================================================
 * Element State Determination
 *
 * Input:  -- MATRIX *d_displ : matrix of global displacements.
 * Output: -- void.
 * =======================================================================
 */

#ifdef  __STDC__
void Elmt_State_Det( MATRIX *d_displ )
#else
void Elmt_State_Det( d_displ )
MATRIX *d_displ;
#endif
{
MATRIX  *dp;
ELEMENT  *el;
ELEMENT_ATTR  *eap;
int  dof_per_elmt;
int  node_no, elmt_no;
int  i, j, k, ii, jj; 
HISTORY_DATA  *hp;
int  UNITS_SWITCH;

#ifdef DEBUG
       printf("*** Enter Elmt_State_Det()\n");
#endif

    UNITS_SWITCH = CheckUnits();

    for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) { 
       el  = &frame->element[elmt_no-1];
       eap = &(frame->eattr[el->elmt_attr_no-1]);

       /* If it is a fiber element, start element state determation */
       /* This block of code should be extended later               */

       if( !(strcmp(eap->elmt_type, "FIBER_2D"))  ||
           !(strcmp(eap->elmt_type, "FIBER_3D"))  ||
           !(strcmp(eap->elmt_type, "FIBER_2DS")) ||
           !(strcmp(eap->elmt_type, "FIBER_3DS")) ) {

          /* Assign element propty */

          for (i = 0; i < NO_ELEMENTS_IN_LIBRARY; i++) { 
             if(!strcmp(elmt_library[i].name, eap->elmt_type)) { 
                k = i;
                break;
             }
          }

          if(array->elmt_no != elmt_no) { /* check if already values assigned */
             array->elmt_no        =  elmt_no;
             array->elmt_attr_no   =  el->elmt_attr_no;
             array->dof_per_node   =  elmt_library[k].no_dof;
             array->no_dimen       =  elmt_library[k].no_dimen;
             array->elmt_type      =  elmt_library[k].name;
             array->nodes_per_elmt = (int) MIN(frame->no_nodes_per_elmt,
                                           elmt_library[k].no_node_per_elmt);
             array->size_of_stiff = array->nodes_per_elmt* array->dof_per_node;

             free((char *) array->material_name);
             array->material_name = SaveString( eap->material );
          }
          array = Assign_p_Array(frame, elmt_no, array, PROPTY);

          /* Assign element coordinates */

          for(j = 1; j <= array->nodes_per_elmt; j++) {
             node_no = el->node_connect[j-1];
             if(node_no != 0) {
                for(i=1;i <= array->no_dimen; i++) {
                   array->coord[i-1][j-1].value = frame->node[node_no -1].coord[i-1].value;
                   if(UNITS_SWITCH == ON )
                      UnitsCopy( array->coord[i-1][j-1].dimen,
                                 frame->node[node_no -1].coord[i-1].dimen );
                }
             }
          }

          /* Transfer fixed displacements */

          for(i=1; i <= array->nodes_per_elmt; i++) {
             k = 1;
             node_no = el->node_connect[i-1];
             for(j = 1; j <= frame->no_dof; j++) {
                jj = frame->node[node_no - 1].bound_id[j-1];  
                if(jj > 0) {
                   array->displ->uMatrix.daa[j-1][i-1] = d_displ->uMatrix.daa[jj-1][0];
                   if( UNITS_SWITCH == ON ) {
                      UnitsCopy(&(array->displ->spRowUnits[j-1]),
                                &(d_displ->spRowUnits[jj-1]));
                      UnitsCopy(&(array->displ->spColUnits[i-1]),
                                &(d_displ->spColUnits[0]));
                   }
                } else {
                   array->displ->uMatrix.daa[j-1][i-1] = frame->node[node_no -1].disp[j-1].value;
                   if( UNITS_SWITCH == ON ) {
                      UnitsCopy(&(array->displ->spRowUnits[j-1]),
                                frame->node[node_no -1].disp[j-1].dimen);
                      UnitsCopy(&(array->displ->spColUnits[i-1]),
                                &(d_displ->spColUnits[0]));
                   }
                }
             }
          }

          /* Retrieve Previous History From The RespondBuffer */

          if( FiberRespondBuffer == (FIBER_RESPONSE *)NULL )
              FatalError("Must have fiber element to use ElmtStateDet()!",
                        (char *) NULL);

          for( ii=0 ; ii < FiberRespondBuffer->total_fiber_elmt ; ++ii ) {
               if( elmt_no == FiberRespondBuffer->history[ii].elmt_no )
                   break;
          }

          hp = &FiberRespondBuffer->history[ii];

          if( !(strcmp(eap->elmt_type, "FIBER_2D")) ||
              !(strcmp(eap->elmt_type, "FIBER_2DS")) )
             Fiber_Elmt_State_Det_2d( array, hp, FiberLoadFlag[ii] );
          else
             Fiber_Elmt_State_Det_3d( array, hp, FiberLoadFlag[ii] );

          el->esp->state = array->elmt_state;
       }
    }

#ifdef DEBUG
       printf("*** Leave Elmt_State_Det()\n");
#endif

}


/*
 * =======================================================================
 * UpdateResponse() : Save element response to frame data structures.
 *
 * Input:  -- void.
 * Output: -- void.
 * =======================================================================
 */

void UpdateResponse() {
ELEMENT             *ep;
ELEMENT_ATTR       *eap;
int             elmt_no;
int        elmt_attr_no;
int                i, j;
int    total_shell_elmt;

SYMBOL             *slp;
int    iInPlaneIntegPts;
int  iThicknessIntegPts;
int       iNO_INTEG_pts;

#ifdef DEBUG
       printf("*** Enter UpdateResponse() \n");
#endif

  total_shell_elmt = 0;

  for(i=1; i<=frame->no_elements; i++) {

    elmt_attr_no = frame->element[i-1].elmt_attr_no;
    eap = &(frame->eattr[elmt_attr_no-1]);

    /*
     * -------------------------------------------------------------------
     * Fiber elements exist, update the load history.
     * Including the flags, stresses and strains for each fiber element.
     * The updated load history will be assign to frame.
     * The related details are in elmt_fiber.c and be done in elmt_fiber.c
     * -------------------------------------------------------------------
     */

    if( !(strcmp(eap->elmt_type, "FIBER_2D"))  ||
        !(strcmp(eap->elmt_type, "FIBER_3D"))  ||
        !(strcmp(eap->elmt_type, "FIBER_2DS")) ||
        !(strcmp(eap->elmt_type, "FIBER_3DS")) ) {
          ep = &(frame->element[i-1]);
          j  = eap->work_fiber->no_fiber +
               (frame->no_dimen-1)*eap->work_fiber->no_shear;
          SaveFiberRespondBuffer( ep, i, frame->no_integ_pt+2, j );
    }

    if( !(strcmp(eap->elmt_type, "SHELL_4N")) ||
        !(strcmp(eap->elmt_type, "SHELL_8N")) )
        total_shell_elmt++;
  }

  /* Update load history for SHELL_4N or SHELL_8N elements */

  if( total_shell_elmt != 0 ) {

      /* no of integration pts in plane/surface */

      slp = lookup("InPlaneIntegPts");  
      if(slp == NULL)
         iInPlaneIntegPts = 2*2;        /* 2x2 as default */
      else
         iInPlaneIntegPts = (int) slp->u.q->value;

      /* no integration pts in thickness direction */

      slp = lookup("ThicknessIntegPts"); 
      if(slp == NULL)
         iThicknessIntegPts = 2;        /* 2 as default */
      else
         iThicknessIntegPts = (int) slp->u.q->value;

      iNO_INTEG_pts = iInPlaneIntegPts*iThicknessIntegPts;

      for(elmt_no = 1; elmt_no <= frame->no_elements; elmt_no++) {
          ep = &frame->element[elmt_no-1];
          save_action(array, ep, &frame->eattr[ep->elmt_attr_no -1],
                      elmt_no, iNO_INTEG_pts);
      }
  }

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

}


/*
 * =======================================================================
 * Save Action ()
 *
 * Input:  -- ARRAY    *p       --
 *         -- ELEMENT *ep       --
 *         -- ELEMENT_ATTR *eap --
 *         -- int elmt_no       --
 *         -- int iNO_INTEG_pts --
 * Output: -- void.
 * =======================================================================
 */

#ifdef __STDC__
void save_action(ARRAY *p, ELEMENT *ep, ELEMENT_ATTR *eap,
                 int elmt_no, int iNO_INTEG_pts )
#else
void save_action(p, ep, eap, elmt_no, iNO_INTEG_pts)
ELEMENT          *ep;
ELEMENT_ATTR    *eap;
ARRAY             *p;
int          elmt_no;
int    iNO_INTEG_pts;

#endif
{
char    *name;
int   i, j, k;
DIMENSIONS *d;

#ifdef DEBUG
       printf("*** Enter save_actions() \n");
#endif

  name = eap->elmt_type;  
  for(i = 1; i <= frame->no_dof; i++ ) {
  for(j = 1; j <= frame->no_nodes_per_elmt; j++) {
       k = frame->no_dof*(j-1) + i;
       ep->rp->Forces->uMatrix.daa[i-1][j-1]
       = RespondBuff[elmt_no-1].Forces->uMatrix.daa[i-1][j-1];
       ep->rp->displ->uMatrix.daa[i-1][j-1] 
       = RespondBuff[elmt_no-1].displ->uMatrix.daa[i-1][j-1];
  }
  }

  ep->rp->max_moment.value = RespondBuff[elmt_no-1].max_moment.value;
  ep->esp->state           = ElmtStateBuff[elmt_no-1];

  /* state, strains parameters */

  switch(ep->esp->state) {
     case 0:   /* ELASTIC state    */
          for(j = 1; j <= iNO_INTEG_pts; j++) {
          ep->rp->effect_pl_strain[j-1] = 0.0;
          for(i = 1; i <= 9; i++) 
              ep->rp->stress->uMatrix.daa[i-1][j-1] 
              = RespondBuff[elmt_no-1].stress->uMatrix.daa[i-1][j-1];
          }
          break;
    case 1:   /* Perfectly plastic or  */
              /* elastic plastic state */
         if(LoadCurve[elmt_no-1].name != NULL) 
            ep->LC_ptr->name = SaveString(LoadCurve[elmt_no-1].name);
         else 
            ep->LC_ptr->name = (char *)NULL;

         for(j = 1; j <= iNO_INTEG_pts; j++) {
            ep->rp->effect_pl_strain[j-1]
               += RespondBuff[elmt_no-1].eff_pl_strain_incr[j-1];

            for(i = 1; i <= 9; i++) {
                ep->rp->stress->uMatrix.daa[i-1][j-1] 
                = RespondBuff[elmt_no-1].stress->uMatrix.daa[i-1][j-1];
                ep->rp->strain_pl->uMatrix.daa[i-1][j-1] 
                += RespondBuff[elmt_no-1].strain_pl_incr->uMatrix.daa[i-1][j-1];
            }
            ep->LC_ptr->R[j-1] = LoadCurve[elmt_no-1].R[j-1];
            ep->LC_ptr->H[j-1] = LoadCurve[elmt_no-1].H[j-1];

            for(i = 1; i <= 6; i++) 
                ep->LC_ptr->back_stress[i-1][j-1] =
                    LoadCurve[elmt_no-1].back_stress[i-1][j-1];
         }
         break;
    default:
         break;
  }

  if(CheckUnits() == ON) {
       switch(CheckUnitsType()) {
         case SI:
              d = DefaultUnits("Pa");
              break;
         case US:
              d = DefaultUnits("psi");
              break;
       }
       for(i = 1; i <= 9; i++) 
           UnitsCopy( &(ep->rp->stress->spRowUnits[i-1]), d );
       free((char *) d->units_name);
       free((char *) d);
  }

#ifdef DEBUG
        printf("******Leaving save_action(): \n");
#endif

}


/*
 *  ===================================================
 *  SetUpRespondBuffer() : Create response buffer......
 *
 *  Input:  -- void.
 *  Output: -- void.
 *  ===================================================
 */

#ifdef __STDC__
void SetUpRespondBuffer()