elmt_fiber2d.c

有限元程序

#else
MATRIX *Rigid_Body_Rotation_2d( length )
QUANTITY length;
#endif
{
int  i, j;
MATRIX *R;

    R  = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 3, 6 );

    for( i=0 ; i < R->iNoRows ; ++i ) {
       for( j=0 ; j < R->iNoColumns ; ++j ) {
          R->uMatrix.daa[i][j] = 0.0;
       }
    }
    R->uMatrix.daa[0][0] = -1.0;
    R->uMatrix.daa[0][3] =  1.0;
    R->uMatrix.daa[1][1] =  1.0/length.value;
    R->uMatrix.daa[1][4] = -1.0/length.value;
    R->uMatrix.daa[1][2] =  1.0;
    R->uMatrix.daa[2][1] =  1.0/length.value;
    R->uMatrix.daa[2][4] = -1.0/length.value;
    R->uMatrix.daa[2][5] =  1.0;

    return( R );
}


/*
 *  =================================================================
 *  Element Transformation (2d case)
 *
 *  Input  : QUANTITY **coord  --
 *         : QUANTITY length   --
 *  Output : MATRIX            --
 *  =================================================================
 */

#ifdef  __STDC__
MATRIX *Element_Transformation_2d( QUANTITY **coord, QUANTITY length )
#else
MATRIX *Element_Transformation_2d( coord, length )
QUANTITY **coord;
QUANTITY length;
#endif
{
MATRIX *Lele;
double cs, sn;
double temp;
int    i, j, k;

    /* Lele3x6 = Rigid3x6 * Trans6x6 */

    cs = (coord[0][1].value - coord[0][0].value)/length.value;
    sn = (coord[1][1].value - coord[1][0].value)/length.value;
    Lele = Rigid_Body_Rotation_2d( length );

    if( cs == 1.0 )
       return( Lele );

    for( i=0 ; i < 3 ; ++i ) {    /* [Q]3x1 */
       for( j=0 ; j < 6 ; j=j+3 ) {
          k = j + 1;
          temp = cs*Lele->uMatrix.daa[i][j] - sn*Lele->uMatrix.daa[i][k];
          Lele->uMatrix.daa[i][k] = sn*Lele->uMatrix.daa[i][j] + cs*Lele->uMatrix.daa[i][k];
          Lele->uMatrix.daa[i][j] = temp;
       }
    }

    return ( Lele );
}


/* 
 *  =======================================================================================
 *  Stress_Strain_Relationship() : Find fiber stress and tangent values according to strain
 *
 *  Input  : HISTORY_DATA  *hp --
 *         : ARRAY      *array --
 *         : FIBER_ATTR *fiber --
 *         : MATRIX        *ex --
 *         : MATRIX       *dex --
 *         : int no_fiber      --
 *         : int no_shear      --
 *         : int sec           --
 *         : int indx          --
 *  Output : void
 *  =======================================================================================
 */ 

#ifdef  __STDC__
void
Stress_Strain_Relationship( HISTORY_DATA *hp, ARRAY *array, FIBER_ATTR *fiber,
                            MATRIX *ex, MATRIX *dex, int no_fiber,
                            int no_shear, int sec, int indx )
#else
void Stress_Strain_Relationship( hp, array, fiber, ex, dex, no_fiber, no_shear, sec, indx )
HISTORY_DATA *hp;
ARRAY *array;
FIBER_ATTR *fiber;
MATRIX  *ex, *dex;
int  no_fiber, no_shear;
int sec, indx;
#endif
{
int  ifib, i, j;
int  total_fiber;
double  fy, ey, ks, kt;
double  k_x, s_x, e_x;

    total_fiber = no_fiber + (array->no_dimen-1)*no_shear;
    for( ifib=0 ; ifib < total_fiber ; ++ifib ) {
       hp->yielding[sec][ifib]  = array->yielding_saved[sec][ifib];
       hp->pre_range[sec][ifib] = array->pre_range_saved[sec][ifib];
       hp->pre_load[sec][ifib]  = array->pre_load_saved[sec][ifib];

       if( indx == 0 ) {
          if( dex->uMatrix.daa[ifib][0] > 0.0 )
             hp->loading[sec][ifib] = 1;
          else if( dex->uMatrix.daa[ifib][0] < 0.0 )
             hp->loading[sec][ifib] = -1;
          else
             hp->loading[sec][ifib] = hp->pre_load[sec][ifib];
       }

       hp->sr->uMatrix.daa[sec][ifib] = array->sr_saved->uMatrix.daa[sec][ifib];
       hp->er->uMatrix.daa[sec][ifib] = array->er_saved->uMatrix.daa[sec][ifib];
       hp->s0->uMatrix.daa[sec][ifib] = array->s0_saved->uMatrix.daa[sec][ifib];
       hp->e0->uMatrix.daa[sec][ifib] = array->e0_saved->uMatrix.daa[sec][ifib];
    }

    for( ifib=0 ; ifib < total_fiber ; ++ifib ) {
       if( ifib < no_fiber )
       {
          ks = fiber[ifib].Es.value;
          kt = fiber[ifib].Et.value;
          fy = fiber[ifib].fy.value;
          ey = fy/ks;
       }
       else if( ifib < 2*no_fiber )
       {
          if( no_shear==1 ) {
             ks = fiber[ifib].Es.value;
             kt = fiber[ifib].Et.value;
             fy = fiber[ifib].fy.value;
             ey = fy/ks;
          }
          else {
             ks = fiber[ifib-no_fiber].Gs.value;
             kt = fiber[ifib-no_fiber].Gt.value;
             fy = fiber[ifib-no_fiber].fv.value;
             ey = fy/ks;
          }
       }
       else
       {
          ks = fiber[ifib-no_fiber-no_shear].Gs.value;
          kt = fiber[ifib-no_fiber-no_shear].Gt.value;
          fy = fiber[ifib-no_fiber-no_shear].fv.value;
          ey = fy/ks;
       }

       if( hp->yielding[sec][ifib] == 0 )  /* first elastic or plastic */
       {
          if( ABS(ex->uMatrix.daa[sec][ifib]) <= ey ) /* first elastic */
          {
             k_x = ks;
             s_x = 0.0;
             e_x = 0.0;
             hp->yielding[sec][ifib] = 0;
             hp->pre_range[sec][ifib] = 0;
          } /* end of first elastic */
          else  /* first from elastic -> plastic, first start yielding */
          {
             k_x = kt;
             hp->s0->uMatrix.daa[sec][ifib] = fy*hp->loading[sec][ifib];
             hp->e0->uMatrix.daa[sec][ifib] = ey*hp->loading[sec][ifib];
             s_x = hp->s0->uMatrix.daa[sec][ifib];
             e_x = hp->e0->uMatrix.daa[sec][ifib];
             hp->yielding[sec][ifib] = 1;
             hp->pre_range[sec][ifib] = 1;
             array->elmt_state = 1;
          } /* end of first start yielding */
       } /* end of yielding[sec][ifib]==0, first elastic or plastic */

       else  /* yielding==1, plastic residual will occurs */
       {
          if( hp->pre_load[sec][ifib] != hp->loading[sec][ifib] )  /* load reversed */
          {
             if( hp->pre_range[sec][ifib] == 1 )
             {
                hp->sr->uMatrix.daa[sec][ifib] = array->sx_saved->uMatrix.daa[sec][ifib];
                hp->er->uMatrix.daa[sec][ifib] = array->ex_saved->uMatrix.daa[sec][ifib];
                k_x = ks;
                s_x = hp->sr->uMatrix.daa[sec][ifib];
                e_x = hp->er->uMatrix.daa[sec][ifib];
                hp->pre_range[sec][ifib] = 0;
             } /* end of pre_range==1 */
             else  /* pre_range==0 */
             {
                if( hp->loading[sec][ifib]*ex->uMatrix.daa[sec][ifib]
                <=  hp->loading[sec][ifib]*hp->er->uMatrix.daa[sec][ifib] )
                {
                   k_x = ks;
                   s_x = hp->sr->uMatrix.daa[sec][ifib];
                   e_x = hp->er->uMatrix.daa[sec][ifib];
                   hp->pre_range[sec][ifib] = 0;
                }
                else
                {
                   if( hp->loading[sec][ifib]*array->ex_saved->uMatrix.daa[sec][ifib]
                   >=  hp->loading[sec][ifib]*hp->er->uMatrix.daa[sec][ifib] )
                   {
                      k_x = ks;
                      s_x = hp->sr->uMatrix.daa[sec][ifib];
                      e_x = hp->er->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 0;
                   }
                   else
                   {
                      k_x = kt;
                      s_x = hp->sr->uMatrix.daa[sec][ifib];
                      e_x = hp->er->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 1;
                   }
                }
             } /* end of pre_range==0 */
          }  /* end of pre_load != loading, load reversed */

          else  /* pre_load==loading, add load in same direction */
          {
             if( hp->pre_range[sec][ifib] == 1 )
             {
                k_x = kt;
                s_x = hp->s0->uMatrix.daa[sec][ifib];
                e_x = hp->e0->uMatrix.daa[sec][ifib];
                hp->pre_range[sec][ifib] = 1;
             }
             else  /* pre_range=0 */
             {
                if( hp->loading[sec][ifib]*array->ex_saved->uMatrix.daa[sec][ifib]
                <=  hp->loading[sec][ifib]*hp->er->uMatrix.daa[sec][ifib] )
                {
                   if( hp->loading[sec][ifib]*ex->uMatrix.daa[sec][ifib]
                   <=  hp->loading[sec][ifib]*hp->er->uMatrix.daa[sec][ifib] )
                   {
                      k_x = ks;
                      s_x = hp->sr->uMatrix.daa[sec][ifib];
                      e_x = hp->er->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 0;
                   }
                   else
                   {
                      k_x = kt;
                      s_x = hp->sr->uMatrix.daa[sec][ifib];
                      e_x = hp->er->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 1;
                   }
                }
                else
                {
                   if( ABS(ex->uMatrix.daa[sec][ifib] - hp->er->uMatrix.daa[sec][ifib]) <= 2*ey )
                   {
                      k_x = ks;
                      s_x = hp->sr->uMatrix.daa[sec][ifib];
                      e_x = hp->er->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 0;
                   }
                   else
                   {
                      hp->s0->uMatrix.daa[sec][ifib] =
                         hp->sr->uMatrix.daa[sec][ifib] + hp->loading[sec][ifib]*2*fy;
                      hp->e0->uMatrix.daa[sec][ifib] =
                         hp->er->uMatrix.daa[sec][ifib] + hp->loading[sec][ifib]*2*ey;
                      k_x = kt;
                      s_x = hp->s0->uMatrix.daa[sec][ifib];
                      e_x = hp->e0->uMatrix.daa[sec][ifib];
                      hp->pre_range[sec][ifib] = 1;
                   }
                }
             } /* end of pre_range=0 */
          }  /* end of pre_load == loading */

       } /* end of yielding==1 */

       hp->tangent->uMatrix.daa[sec][ifib] = k_x;
       hp->stress->uMatrix.daa[sec][ifib]  = s_x + k_x*(ex->uMatrix.daa[sec][ifib]-e_x);
    }
}


/*
 *  =========================================================
 *  Get Gauss-Lobatto abscissas and weights
 * 
 *  Input : abscissas = section position in one element
 *        : weights   = array of weighting coefficients.
 *        : n = integration points = number of section in one
 *              element
 *  Output : void.
 *  =========================================================
 */

#ifdef  __STDC__
void Gauss_Lobatto( double *abscissas, double *weights , int n )
#else
void Gauss_Lobatto( abscissas, weights , n )
double *abscissas, *weights;
int n;
#endif
{
	/* Gauss-Lobatto integration */
	/* integral(-1,1)f(x)dx = A*f(-1) + sum(1,n)Ai*f(xi) + A*f(1) */

	switch(n)
	{
	   case 2:
	      abscissas[0] = -1.0;
	      abscissas[1] = -0.4472135954;
	      abscissas[2] =  0.4472135954;
	      abscissas[3] =  1.0;

	      weights[0] = 0.1666666666;
	      weights[1] = 0.8333333333;
	      weights[2] = 0.8333333333;
	      weights[3] = 0.1666666666;
	      break;

	   case 3:
	      abscissas[0] = -1.0;
	      abscissas[1] = -0.6546536707;
	      abscissas[2] =  0.0;
	      abscissas[3] =  0.6546536707;
	      abscissas[4] =  1.0;

	      weights[0] = 0.1000000000;
	      weights[1] = 0.5444444444;
	      weights[2] = 0.7111111111;
	      weights[3] = 0.5444444444;
	      weights[4] = 0.1000000000;
	      break;

	   case 4:
	      abscissas[0] = -1.0;
	      abscissas[1] = -0.7650553239;
	      abscissas[2] = -0.2852315164;
	      abscissas[3] =  0.2852315164;
	      abscissas[4] =  0.7650553239;
	      abscissas[5] =  1.0;

	      weights[0] = 0.0666666666;
	      weights[1] = 0.3784749562;
	      weights[2] = 0.5548583770;
	      weights[3] = 0.5548583770;
	      weights[4] = 0.3784749562;
	      weights[5] = 0.0666666666;
	      break;

	   case 5:
	      abscissas[0] = -1.0;
	      abscissas[1] = -0.8302238962;
	      abscissas[2] = -0.4688487934;
	      abscissas[3] =  0.0;
	      abscissas[4] =  0.4688487934;
	      abscissas[5] =