elmt_fiber2d.c

有限元程序

#endif

    /* Initial Residual Deformation Matrix, Resisting Force Matrix */

    s   = MatrixAllocIndirect(   "s",  DOUBLE_ARRAY, 3, 1 );
    DRx = MatrixAllocIndirect( "DRx",  DOUBLE_ARRAY, 3, 1 );

    /* Element Converge, j ; At Least Do Once */

    do
    {
       /* Reset Flexibility Matrix F, And Residual Matrix s */

       for( ii=0 ; ii < F->iNoRows ; ++ii )
       for( jj=0 ; jj < F->iNoColumns ; ++jj )
          F->uMatrix.daa[ii][jj] = 0.0;

       for( ii=0 ; ii < s->iNoRows ; ++ii ) {
          s->uMatrix.daa[ii][0] = 0.0;
       }

       /* Section State Determination */

       dQ = MatrixMult( K, dq );
       MatrixAddReplace( Q, dQ );

#ifdef DEBUG2
printf("\n K in the do loop\n");
MatrixPrintVar( K, (MATRIX *)NULL );
printf("\n dq\n");
MatrixPrintVar( dq, (MATRIX *)NULL );
printf("\n dQ\n");
MatrixPrintVar( dQ, (MATRIX *)NULL );
#endif

       for( sec=0 ; sec < no_section; ++sec ) {
          xi = saved_data[sec].xi;
          scale = saved_data[sec].wi;
          fx = saved_data[sec].fx;
          Dx = saved_data[sec].Dx;
          rx = saved_data[sec].rx;

#ifdef DEBUG2
printf("\n sec = %i\n",sec); 
printf("\n fx at the beginnig of do loop\n");
MatrixPrintVar( fx, (MATRIX *)NULL );
printf("\n Dx\n");
MatrixPrintVar( Dx, (MATRIX *)NULL );
#endif

          Force_Interpolation_Matrix_2d( bx, xi, length );
          Linear_Geometric_Matrix_2d( lx, fiber, no_fiber, no_shear );

          dDx = MatrixMult( bx, dQ );
          MatrixAddReplace( Dx, dDx );

          temp_m1 = MatrixMult( fx, dDx );
          ddx = MatrixAdd( rx, temp_m1 );
          dex = MatrixMult( lx, ddx );
          MatrixFree( temp_m1 );
          for( ifib=0 ; ifib < total_fiber ; ++ifib )
               ex->uMatrix.daa[sec][ifib] += dex->uMatrix.daa[ifib][0];

#ifdef DEBUG2
printf("\n dDx\n");
MatrixPrintVar( dDx, (MATRIX *)NULL );
printf("\n ddx\n");
MatrixPrintVar( ddx, (MATRIX *)NULL );
printf("\n dex\n");
MatrixPrintVar( dex, (MATRIX *)NULL );
#endif

          Stress_Strain_Relationship( hp, p, fiber, ex, dex, no_fiber, no_shear, sec, flag[sec] );
          flag[sec]++;

#ifdef DEBUG2
printf("\n stress\n");
MatrixPrintVar( stress, (MATRIX *)NULL );
printf("\n tangent\n");
MatrixPrintVar( tangent, (MATRIX *)NULL );
#endif

          MatrixFree( dDx );
          MatrixFree( ddx );
          MatrixFree( dex );

          Section_Tangent_Stiffness_2d( kx, fiber, no_fiber, no_shear, sec, elmt_no );
          MatrixFree( fx );
          fx = MatrixInverse( kx );
          saved_data[sec].fx = fx;

          Section_Resisting_Force_2d( DRx, stress, fiber, no_fiber, no_shear, sec );
          DUx = MatrixSub( Dx, DRx );
          MatrixFree( rx );
          rx = MatrixMult( fx, DUx );
          saved_data[sec].rx = rx;
          MatrixFree( DUx );

#ifdef DEBUG2
printf("\n fx after Tangent\n");
MatrixPrintVar( fx, (MATRIX *)NULL );
printf("\n rx\n");
MatrixPrintVar( rx, (MATRIX *)NULL );
#endif

          bxtrans = MatrixTranspose( bx );
          temp_m1 = MatrixMult( bxtrans, fx );
          temp_m2 = MatrixMult( temp_m1, bx );
          MatrixFree( temp_m1 );
          temp_m1 = MatrixScale( temp_m2, scale );
          MatrixAddReplace( F, temp_m1 );
          MatrixFree( temp_m1 );
          MatrixFree( temp_m2 );

          temp_m1 = MatrixMult( bxtrans, rx );
          temp_m2 = MatrixScale( temp_m1, scale );
          MatrixAddReplace( s, temp_m2 );
          MatrixFree( temp_m1 );
          MatrixFree( temp_m2 );
          MatrixFree( bxtrans );
       }
       MatrixFree( K );
       K = MatrixInverse( F );

       for( ii=0 ; ii < dq->iNoRows ; ++ii )
          dq->uMatrix.daa[ii][0] = -s->uMatrix.daa[ii][0];

       MatrixFree( dQ );

    }  while ( dVmatrixL2Norm(s->uMatrix.daa, s->iNoRows, s->iNoColumns) > 0.001 );

    for( sec=0 ; sec < no_section ; ++sec ) {
       MatrixFree( saved_data[sec].fx );
       MatrixFree( saved_data[sec].Dx );
       MatrixFree( saved_data[sec].rx );
    }
    free((char *) saved_data);
    free((char *) abscissas);
    free((char *) weights);
    MatrixFree( dpe );
    MatrixFree( kx );
    MatrixFree( bx );
    MatrixFree( lx );
    MatrixFree( F );
    MatrixFree( K );
    MatrixFree( L );
    MatrixFree( s );
    MatrixFree( dq );
    MatrixFree( DRx );

    if( UNITS_SWITCH == ON )
       SetUnitsOn();

#ifdef DEBUG
    printf("*** ===============================\n");
    printf("*** LEAVE Fiber_Elmt_State_Det_2d()\n");
    printf("*** ===============================\n");
#endif

}



/*
 *  =================================================================
 *  Force Interpolation Matrix (2d case)
 *
 *  Input  : Array  *bx -- pointer to .....
 *         : double x   -- ....
 *         : QUANTITY L -- 
 *  Output : void
 *  =================================================================
 */



#ifdef  __STDC__
void Force_Interpolation_Matrix_2d( MATRIX *bx, double x, QUANTITY L )
#else
void Force_Interpolation_Matrix_2d( bx, x, L )
MATRIX *bx;
double   x;
QUANTITY L;
#endif
{
    bx->uMatrix.daa[0][0] = 1.0;
    bx->uMatrix.daa[0][1] = 0.0;
    bx->uMatrix.daa[0][2] = 0.0;
    bx->uMatrix.daa[1][0] = 0.0;
    bx->uMatrix.daa[1][1] = x/L.value-1.0;
    bx->uMatrix.daa[1][2] = x/L.value;
    bx->uMatrix.daa[2][0] = 0.0;
    bx->uMatrix.daa[2][1] = 1.0/L.value;
    bx->uMatrix.daa[2][2] = 1.0/L.value;
}


/*
 *  =================================================================
 *  Linear Geometric Matrix (2d case).
 *
 *  Input  : MATRIX *lx        --
 *         : FIBER_ATTR *fiber --
 *         : int no_fiber      --
 *           int no_shear      --
 *  Output : void
 *  =================================================================
 */

#ifdef  __STDC__
void Linear_Geometric_Matrix_2d( MATRIX *lx, FIBER_ATTR *fiber, int no_fiber, int no_shear )
#else
void Linear_Geometric_Matrix_2d( lx, fiber, no_fiber, no_shear )
MATRIX	*lx;
FIBER_ATTR  *fiber;
int  no_fiber, no_shear;
#endif
{
int	ifib, total_fiber;

    total_fiber = no_fiber + no_shear;
    for( ifib=0 ; ifib < no_fiber ; ++ifib )
    {
       lx->uMatrix.daa[ifib][0] =  1.0;
       lx->uMatrix.daa[ifib][1] = -fiber[ifib].y.value;
       lx->uMatrix.daa[ifib][2] =  0.0;
    }
    for( ifib=no_fiber ; ifib < total_fiber ; ++ifib )
    {
       lx->uMatrix.daa[ifib][0] =  0.0;
       lx->uMatrix.daa[ifib][1] =  0.0;
       lx->uMatrix.daa[ifib][2] =  1.0;
    }
}


/*
 *  =================================================================
 *  Section Tangent Stiffness Matrix
 *
 *  Input  : MATRIX *kx        --
 *         : FIBER_ATTR *fiber --
 *         : int no_fiber      --
 *         : int no_shear      --
 *         : int sec           -- 
 *         : int elmt_no       --
 *  Output : void
 *  =================================================================
 */


#ifdef  __STDC__
void
Section_Tangent_Stiffness_2d( MATRIX *kx, FIBER_ATTR *fiber, int no_fiber,
                              int no_shear, int sec, int elmt_no )
#else
void
Section_Tangent_Stiffness_2d( kx, fiber, no_fiber, no_shear, sec, elmt_no )
MATRIX *kx;
FIBER_ATTR *fiber;
int no_fiber, no_shear, sec;
int elmt_no;
#endif
{
int	  ifib, i, j;
int      total_fiber;
double	y, z, A, Exi;
HISTORY_DATA     *hp;
MATRIX      *tangent;

    /* Zero tangent stiffness matrix */

    for( i=0 ; i < kx->iNoRows ; ++i )
       for( j=0 ; j < kx->iNoColumns ; ++j )
          kx->uMatrix.daa[i][j] = 0.0;

    /* For strain at xi, get E(xi) for the fiber */

    hp = FiberElmtHistory( elmt_no );
    tangent = hp->tangent;

    for( ifib=0 ; ifib < no_fiber ; ++ifib ) {

       y   = fiber[ifib].y.value;
       A   = fiber[ifib].area.value;
       Exi = tangent->uMatrix.daa[sec][ifib];

       kx->uMatrix.daa[0][0] = kx->uMatrix.daa[0][0] + Exi*A;
       kx->uMatrix.daa[0][1] = kx->uMatrix.daa[0][1] - Exi*A*y;
       kx->uMatrix.daa[1][0] = kx->uMatrix.daa[0][1];
       kx->uMatrix.daa[1][1] = kx->uMatrix.daa[1][1] + Exi*A*y*y;
    }

    if( no_shear == 1 ) {
       kx->uMatrix.daa[2][2] = tangent->uMatrix.daa[sec][no_fiber]*fiber[no_fiber].area.value;
    } else {
       for( ifib=no_fiber ; ifib < total_fiber ; ++ifib )
          kx->uMatrix.daa[2][2] = kx->uMatrix.daa[2][2] +
                                  tangent->uMatrix.daa[sec][ifib]*fiber[ifib-no_fiber].area.value;
    }
}


/*
 *  =================================================================
 *  Section Resisting Force
 *
 *  Input  : MATRIX *DRx       --
 *         : MATRIX *stress    --
 *         : FIBER_ATTR *fiber --
 *         : int no_fiber      --
 *         : int no_shear      --
 *         : int sec           --
 *  Output : void
 *  =================================================================
 */

#ifdef  __STDC__
void
Section_Resisting_Force_2d( MATRIX *DRx, MATRIX *stress, FIBER_ATTR *fiber,
                            int no_fiber, int no_shear, int sec )
#else
void
Section_Resisting_Force_2d( DRx, stress, fiber, no_fiber, no_shear, sec )
MATRIX *DRx;
MATRIX *stress;
FIBER_ATTR *fiber;
int  no_fiber, no_shear, sec;
#endif
{
int  ifib, i, j;
int  total_fiber;

    total_fiber = no_fiber + no_shear;
    for( i=0 ; i < DRx->iNoRows ; ++i )
       DRx->uMatrix.daa[i][0] = 0.0;

    for( ifib=0 ; ifib < no_fiber ; ++ifib )
    {
       DRx->uMatrix.daa[0][0] = DRx->uMatrix.daa[0][0] +
            stress->uMatrix.daa[sec][ifib]*fiber[ifib].area.value;
       DRx->uMatrix.daa[1][0] = DRx->uMatrix.daa[1][0] -
            stress->uMatrix.daa[sec][ifib]*fiber[ifib].area.value*fiber[ifib].y.value;
    }

    if( no_shear == 1 ) {
       DRx->uMatrix.daa[2][0] = stress->uMatrix.daa[sec][no_fiber]*fiber[no_fiber].area.value;
    }
    else {
    for( ifib=no_fiber ; ifib < total_fiber ; ++ifib )
       DRx->uMatrix.daa[2][0] = DRx->uMatrix.daa[2][0] +
            stress->uMatrix.daa[sec][ifib]*fiber[ifib-no_fiber].area.value;
    }
}


/*
 *  =================================================================
 *  Rigid Body Rotation
 *
 *  Input  : QUANTITY length  --
 *  Output : MATRIX           --
 *  =================================================================
 */

#ifdef  __STDC__
MATRIX *Rigid_Body_Rotation_2d( QUANTITY length )