fe_nonlinear.c

有限元程序

#else
void SetUpRespondBuffer()
#endif
{
int                i, j;
int         iNoIntegPts; 

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

ELEMENT_ATTR   *eap;
FIBER_ELMT     *fep;

int	total_fiber_elmt;
int	total_shell_elmt;

#ifdef DEBUG
   printf("*** Enter SetUpRespondBuffer() \n");
   printf("*** No_elements = %d           \n", frame->no_elements);
#endif

  total_fiber_elmt = 0;
  total_shell_elmt = 0;

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

     eap = lookup(frame->element[i-1].elmt_attr_name)->u.eap;
     if(eap == NULL)
        FatalError("Elmt_Attribute name not found",(char *)NULL);

     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")) )
           total_fiber_elmt++;

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

  /*
   * ----------------------------------------------------------------------
   * Fiber elements exist, setup the needed space for storing load history.
   * Including the flags, stresses and strains for each fiber element.
   * The storage space will be static and assigned to frame until updated.
   * The related details are in elmt_fiber.c and be done in elmt_fiber.c
   * ----------------------------------------------------------------------
   */

  if( total_fiber_elmt != 0 )
      SetUpFiberRespondBuffer( total_fiber_elmt, frame );

  /* SHELL_4N or SHELL_8N elements exist, setup the needed space for storing load history. */

  if( total_shell_elmt != 0 ) {

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

    /* Number of integration pts in thickness direction*/

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

    iNoIntegPts = iInPlaneIntegPts*iThicknessIntegPts;


    RespondBuff = (RESPONSE *) MyCalloc(frame->no_elements, sizeof(RESPONSE)); 
    ElmtStateBuff = (int *)    MyCalloc(frame->no_elements, sizeof(int)); 
    LoadCurve     = (MATER_LOAD_CURVE *) MyCalloc(frame->no_elements,
                                           sizeof(MATER_LOAD_CURVE)); 

    for(i = 1; i <= frame->no_elements; i++) {
        RespondBuff[i-1].Forces 
          = MatrixAllocIndirect("NodalForce",DOUBLE_ARRAY,6,frame->no_nodes);
        RespondBuff[i-1].displ
          = MatrixAllocIndirect("NodalDispl",DOUBLE_ARRAY,6,frame->no_nodes);
        RespondBuff[i-1].stress
          = MatrixAllocIndirect("Stress",DOUBLE_ARRAY,9,iNoIntegPts);
        RespondBuff[i-1].strain_pl
          = MatrixAllocIndirect("PlasticStrain",DOUBLE_ARRAY,9,iNoIntegPts);
        RespondBuff[i-1].strain_pl_incr
          = MatrixAllocIndirect("PlasticStrainIncr",DOUBLE_ARRAY,9,iNoIntegPts);
        RespondBuff[i-1].effect_pl_strain
          = (double *) MyCalloc(iNoIntegPts, sizeof(double));
        RespondBuff[i-1].eff_pl_strain_incr 
          = (double *) MyCalloc(iNoIntegPts, sizeof(double));

        LoadCurve[i-1].name = (char *) MyCalloc(1,sizeof(char));
        LoadCurve[i-1].R    = (double *) MyCalloc(iNoIntegPts, sizeof(double));
        LoadCurve[i-1].H    = (double *) MyCalloc(iNoIntegPts, sizeof(double));
        LoadCurve[i-1].back_stress = MatrixAllocIndirectDouble(6, iNoIntegPts);
    }
  }  /* end of total_shell_elmt != 0, SHELL_4N or SHELL_8N elements exist */

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

}

/*
 * =======================================================================
 * SaveRespondBuffer() : Save response buffer..... 
 *
 * Input:  -- ARRAy *p     --
 *         -- int integ_pt --
 * Output: -- void.
 * =======================================================================
 */

#ifdef __STDC__
void SaveRespondBuffer(ARRAY *p, int integ_pt) 
#else
void SaveRespondBuffer(p, integ_pt) 
ARRAY     *p;
int integ_pt;
#endif
{
int i, j, k, kk;
      
#ifdef DEBUG
       printf("*** Enter SaveRespondBuffer(): \n");
#endif

    kk = integ_pt;
    ElmtStateBuff[p->elmt_no-1] = p->elmt_state;

    if(p->LC_ptr->name != NULL) 
       LoadCurve[p->elmt_no-1].name = SaveString(p->LC_ptr->name);
     else
       LoadCurve[p->elmt_no-1].name = NULL;
    
    LoadCurve[p->elmt_no-1].R[kk-1] = p->LC_ptr->R[kk-1];
    LoadCurve[p->elmt_no-1].H[kk-1] = p->LC_ptr->H[kk-1];

    for(i = 1; i <= 6; i++)
        LoadCurve[p->elmt_no-1].back_stress[i-1][kk-1]
        = p->LC_ptr->back_stress[i-1][kk-1];

    RespondBuff[p->elmt_no-1].effect_pl_strain[kk-1]
    = p->effect_pl_strain[kk-1];

    RespondBuff[p->elmt_no-1].eff_pl_strain_incr[kk-1]
    = p->eff_pl_strain_incr[kk-1];

    for(i = 1; i <= p->dof_per_node; i++) {
        RespondBuff[p->elmt_no-1].stress->uMatrix.daa[i-1][kk-1]
        = p->stress->uMatrix.daa[i-1][kk-1];
        RespondBuff[p->elmt_no-1].strain_pl->uMatrix.daa[i-1][kk-1]
        = p->strain_pl->uMatrix.daa[i-1][kk-1];
        RespondBuff[p->elmt_no-1].strain_pl_incr->uMatrix.daa[i-1][kk-1]
        = p->strain_pl_incr->uMatrix.daa[i-1][kk-1];
    }

    for(i = 1; i <= p->dof_per_node; i++ ) {
    for(j = 1; j <= p->nodes_per_elmt; j++) {
        k = p->dof_per_node*(j-1) + i;
        RespondBuff[p->elmt_no-1].Forces->uMatrix.daa[i-1][j-1] =
                    p->nodal_loads[k-1].value;
        RespondBuff[p->elmt_no-1].displ->uMatrix.daa[i-1][j-1] =
                    p->displ->uMatrix.daa[i-1][j-1];
    }
    }

    /* What's this ??? */

    if((p->elmt_type != NULL) && !strcmp(p->elmt_type, "FRAME_2D"))
        RespondBuff[p->elmt_no-1].max_moment.value 
           = MAX( p->nodal_loads[2].value, p->nodal_loads[5].value);
    if((p->elmt_type != NULL) && !strcmp(p->elmt_type, "FRAME_3D"))
        RespondBuff[p->elmt_no-1].max_moment.value
           = MAX( p->nodal_loads[5].value, p->nodal_loads[11].value);

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

}

/*
 * =======================================================================
 * SaveFiberRespondBuffer() : Save the static flags, stresses and strains
 *                            of each fiber element to frame.
 *
 * Input:  -- ELEMENT *el,
 *         -- int elmt_no,
 *         -- int no_section,
 *         -- int total_no_fiber
 * Output: -- void.
 * =======================================================================
 */

#ifdef  __STDC__
void SaveFiberRespondBuffer( ELEMENT *el, int elmt_no,
                             int no_section, int total_no_fiber )
#else
void SaveFiberRespondBuffer( el, elmt_no, no_section, total_no_fiber )
ELEMENT  *el;
int  elmt_no, no_section, total_no_fiber;
#endif
{
HISTORY_DATA      *hp;
int     ii, ifib, sec;

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

    hp = &FiberRespondBuffer->history[ii];

    for( sec=0 ; sec < no_section ; ++sec ) {
    for( ifib=0 ; ifib < total_no_fiber ; ++ifib ) {
       hp->pre_load[sec][ifib] = hp->loading[sec][ifib];
       el->esp->yielding_saved[sec][ifib]  = hp->yielding[sec][ifib];
       el->esp->pre_range_saved[sec][ifib] = hp->pre_range[sec][ifib];
       el->esp->pre_load_saved[sec][ifib]  = hp->pre_load[sec][ifib];

       el->rp->sr_saved->uMatrix.daa[sec][ifib] = hp->sr->uMatrix.daa[sec][ifib];
       el->rp->er_saved->uMatrix.daa[sec][ifib] = hp->er->uMatrix.daa[sec][ifib];
       el->rp->s0_saved->uMatrix.daa[sec][ifib] = hp->s0->uMatrix.daa[sec][ifib];
       el->rp->e0_saved->uMatrix.daa[sec][ifib] = hp->e0->uMatrix.daa[sec][ifib];
       el->rp->sx_saved->uMatrix.daa[sec][ifib] = hp->stress->uMatrix.daa[sec][ifib];
       el->rp->ex_saved->uMatrix.daa[sec][ifib] = hp->strain->uMatrix.daa[sec][ifib];
    }
    }

    /* Initial flag matrix, for each load step */
    for( sec=0 ; sec < no_section; ++sec )
       FiberLoadFlag[ii][sec] = 0;
}

/*
 * ========================================================================
 * SetUpFiberRespondBuffer() : Setup the static flags, stresses and strains
 *                             to store load history. 
 *
 * Input:  -- int total_fiber_elmt 
 *         -- EFRAME *frp 
 * Output: -- void.
 * ========================================================================
 */

#ifdef  __STDC__
void SetUpFiberRespondBuffer( int total_fiber_elmt, EFRAME *frp )
#else
void SetUpFiberRespondBuffer( total_fiber_elmt, frp )
int  total_fiber_elmt;
EFRAME *frp;
#endif
{
HISTORY_DATA    *hp;
int      ii, jj, kk;
int       ifib, sec;
int         elmt_no;
int      no_section;
int        no_fiber;
int        no_shear;
int  total_no_fiber;
int    elmt_attr_no;
int    UNITS_SWITCH;
QUANTITY         Es;
DIMENSIONS   *dimen;

    UNITS_SWITCH = CheckUnits();
    if( UNITS_SWITCH == ON )
       SetUnitsOff();

    FiberRespondBuffer = (FIBER_RESPONSE *)MyCalloc(1, sizeof(FIBER_RESPONSE));
    FiberRespondBuffer->total_fiber_elmt = total_fiber_elmt;
    FiberRespondBuffer->history = (HISTORY_DATA *) MyCalloc( total_fiber_elmt,
                                                             sizeof(HISTORY_DATA));

    no_section = frp->no_integ_pt + 2;  /* include 2 end sections */

    /* Initial flag matrix, for each load step */

    FiberLoadFlag = (int **)MyCalloc( total_fiber_elmt, sizeof(int *) );
    for( ii=0 ; ii < total_fiber_elmt; ++ii ) {
       FiberLoadFlag[ii] = (int *)MyCalloc( no_section, sizeof(int) );
       for( sec=0 ; sec < no_section; ++sec )
          FiberLoadFlag[ii][sec] = 0;
    }

    jj = 0;
    for( ii=1 ; ii <= frp->no_elements; ++ii ) {

       elmt_attr_no = frp->element[ii-1].elmt_attr_no;
       if( !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_2D"))
        || !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_3D"))
        || !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_2DS"))
        || !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_3DS")) )
       { 
	  /* total no. of fiber includes shear fibers */
	  no_fiber = frp->eattr[elmt_attr_no-1].work_fiber->no_fiber;
	  no_shear = frp->eattr[elmt_attr_no-1].work_fiber->no_shear;
	  total_no_fiber = no_fiber + (frp->no_dimen-1)*no_shear;
          elmt_no = ii;
          jj++;
          hp = &FiberRespondBuffer->history[jj-1];
          hp->elmt_no = elmt_no;
          hp->sr = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                  no_section, total_no_fiber);
          hp->er = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                  no_section, total_no_fiber);
          hp->s0 = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                  no_section, total_no_fiber);
          hp->e0 = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                  no_section, total_no_fiber);
          hp->stress = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                       no_section, total_no_fiber);
          hp->strain = (MATRIX *) MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                      no_section, total_no_fiber);
          hp->tangent = (MATRIX *)MatrixAllocIndirect((char *)NULL, DOUBLE_ARRAY,
                                                      no_section, total_no_fiber);

          /* Initial the fiber tangent value into kx */

          for( ifib=0 ; ifib < no_fiber ; ++ifib ) {
             Es = frp->eattr[elmt_attr_no-1].work_fiber->fiber[ifib].Es;
             for( sec=0 ; sec < no_section ; ++sec )
                hp->tangent->uMatrix.daa[sec][ifib] = Es.value;
          }

          if( !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_2D")) 
           || !(strcmp(frp->eattr[elmt_attr_no-1].elmt_type, "FIBER_3D")) ) {

             for( ifib=no_fiber ; ifib < total_no_fiber ; ++ifib ) {
                Es = frp->eattr[elmt_attr_no-1].work_fiber->fiber[ifib].Es;
                for( sec=0 ; sec < no_section ; ++sec )
                   hp->tangent->uMatrix.daa[sec][ifib] = Es.value;
             }

          } else {

             for( ifib=no_fiber ; ifib < (no_fiber+no_shear) ; ++ifib ) {
                Es = frp->eattr[elmt_attr_no-1].work_fiber->fiber[ifib-no_fiber].Gs;
                for( sec=0 ; sec < no_section ; ++sec )
                   hp->tangent->uMatrix.daa[sec][ifib] = Es.value;
             }
             for( ifib=(no_fiber+no_shear) ; ifib < total_no_fiber ; ++ifib ) {
                Es = frp->eattr[elmt_attr_no-1].work_fiber->fiber[ifib-no_fiber-no_shear].Gs;
                for( sec=0 ; sec < no_section ; ++sec )
                   hp->tangent->uMatrix.daa[sec][ifib] = Es.value;
             }

          }

          hp->yielding  = (int **)MyCalloc( no_section, sizeof(int *) );
          hp->pre_load  = (int **)MyCalloc( no_section, sizeof(int *) );
          hp->pre_range = (int **)MyCalloc( no_section, sizeof(int *) );
          hp->loading   = (int **)MyCalloc( no_section, sizeof(int *) );

          for( kk=0 ; kk < no_section ; ++kk ) {

	     hp->yielding[kk]  = (int *)MyCalloc( total_no_fiber, sizeof(int) );
	     hp->pre_load[kk]  = (int *)MyCalloc( total_no_fiber, sizeof(int) );
	     hp->pre_range[kk] = (int *)MyCalloc( total_no_fiber, sizeof(int) );
	     hp->loading[kk]   = (int *)MyCalloc( total_no_fiber, sizeof(int) );

          }
       } /* end of if loop for fiber element */
    } /* end of for loop for all element */

    if( UNITS_SWITCH == ON )
       SetUnitsOn();
}

/*
 * =======================================================================
 * FiberElmtHistory() : Fiber element history .....
 *
 * Input:  -- int elmt_no :
 * Output: -- void.
 * =======================================================================
 */

#ifdef  __STDC__
HISTORY_DATA *FiberElmtHistory( int elmt_no )
#else
HISTORY_DATA *FiberElmtHistory( elmt_no )
int elmt_no;
#endif
{
int  ii;
HISTORY_DATA *hp;

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

    hp = &FiberRespondBuffer->history[ii];

    return( hp );
}