elmt_fiber3d.c

有限元程序

          flag[sec]++;
#ifdef DEBUG
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_3d( kx, fiber, no_fiber, no_shear, sec, elmt_no );
          MatrixFree( fx );
          fx = MatrixInverse( kx );
          saved_data[sec].fx = fx;

          Section_Resisting_Force_3d( 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 DEBUG
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  __STDC__
void Force_Interpolation_Matrix_3d( MATRIX *bx, double x, QUANTITY L )
#else
void Force_Interpolation_Matrix_3d( 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[0][3] = 0.0;
    bx->uMatrix.daa[0][4] = 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[1][3] = 0.0;
    bx->uMatrix.daa[1][4] = 0.0;

    bx->uMatrix.daa[2][0] = 0.0;
    bx->uMatrix.daa[2][1] = 0.0;
    bx->uMatrix.daa[2][2] = 0.0;
    bx->uMatrix.daa[2][3] = x/L.value-1.0;
    bx->uMatrix.daa[2][4] = x/L.value;

    bx->uMatrix.daa[3][0] = 0.0;
    bx->uMatrix.daa[3][1] = 0.0;
    bx->uMatrix.daa[3][2] = 0.0;
    bx->uMatrix.daa[3][3] = 1.0/L.value;
    bx->uMatrix.daa[3][4] = 1.0/L.value;

    bx->uMatrix.daa[4][0] = 0.0;
    bx->uMatrix.daa[4][1] = -1.0/L.value;
    bx->uMatrix.daa[4][2] = -1.0/L.value;
    bx->uMatrix.daa[4][3] = 0.0;
    bx->uMatrix.daa[4][4] = 0.0;
}

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

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

#ifdef  __STDC__
void Section_Tangent_Stiffness_3d( MATRIX *kx, FIBER_ATTR *fiber, int no_fiber, int no_shear, int sec, int elmt_no )
#else
void Section_Tangent_Stiffness_3d( 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;
double	y, z, A, Exi;
HISTORY_DATA     *hp;
MATRIX      *tangent;

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

    /* according to 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;
       z = fiber[ifib].z.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*z;
       kx->uMatrix.daa[0][2] = kx->uMatrix.daa[0][2] - Exi*A*y;
       kx->uMatrix.daa[1][1] = kx->uMatrix.daa[1][1] + Exi*A*z*z;
       kx->uMatrix.daa[1][2] = kx->uMatrix.daa[1][2] - Exi*A*y*z;
       kx->uMatrix.daa[2][2] = kx->uMatrix.daa[2][2] + Exi*A*y*y;
    }
    kx->uMatrix.daa[1][0] = kx->uMatrix.daa[0][1];
    kx->uMatrix.daa[2][0] = kx->uMatrix.daa[0][2];
    kx->uMatrix.daa[2][1] = kx->uMatrix.daa[1][2];

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

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

#ifdef  __STDC__
void Section_Resisting_Force_3d( MATRIX *DRx, MATRIX *stress, FIBER_ATTR *fiber, int no_fiber, int no_shear, int sec )
#else
void Section_Resisting_Force_3d( 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;

    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].z.value;
       DRx->uMatrix.daa[2][0] = DRx->uMatrix.daa[2][0] - stress->uMatrix.daa[sec][ifib]*fiber[ifib].area.value*fiber[ifib].y.value;
    }

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

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

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

    R  = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 5, 12 );

    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][6]  =  1.0;
    R->uMatrix.daa[1][2]  = -1.0/length.value;
    R->uMatrix.daa[1][4]  =  1.0;
    R->uMatrix.daa[1][8]  =  1.0/length.value;
    R->uMatrix.daa[2][2]  = -1.0/length.value;
    R->uMatrix.daa[2][8]  =  1.0/length.value;
    R->uMatrix.daa[2][10] =  1.0;
    R->uMatrix.daa[3][1]  =  1.0/length.value;
    R->uMatrix.daa[3][5]  =  1.0;
    R->uMatrix.daa[3][7]  = -1.0/length.value;
    R->uMatrix.daa[4][1]  =  1.0/length.value;
    R->uMatrix.daa[4][7]  = -1.0/length.value;
    R->uMatrix.daa[4][11] =  1.0;

    return( R );
}

#ifdef  __STDC__
MATRIX *Element_Transformation_3d( QUANTITY **coord, QUANTITY length )
#else
MATRIX *Element_Transformation_3d( coord, length )
QUANTITY **coord;
QUANTITY length;
#endif
{
MATRIX *Lele;
double cx, cy, cz, cxz;
double temp1, temp2, temp3;
int    i, j, k;

    /* Lele5x12 = Rigid5x12 * Trans12x12 */

    cx = (coord[0][1].value - coord[0][0].value)/length.value;
    cy = (coord[1][1].value - coord[1][0].value)/length.value;
    cz = (coord[2][1].value - coord[2][0].value)/length.value;
    cxz = sqrt( cx*cx + cz*cz );
    Lele = Rigid_Body_Rotation_3d( length );

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

    /* rotation of vertical axes */
    if( cxz == 0.0 ) {
       for( i=0 ; i < 5 ; ++i ) {  /* [Q]5x1 */
          for( j=0 ; j < 12 ; j=j+3 ) {
             temp1 = Lele->uMatrix.daa[i][j];
             temp2 = Lele->uMatrix.daa[i][j+1];
             Lele->uMatrix.daa[i][j]   = -temp2*cy;
             Lele->uMatrix.daa[i][j+1] =  temp1*cy;
          }
       }
       return ( Lele );
    }

    for( i=0 ; i < 5 ; ++i ) {  /* [Q]5x1 */
       for( j=0 ; j < 12 ; j=j+3 ) {
          temp1 = Lele->uMatrix.daa[i][j];
          temp2 = Lele->uMatrix.daa[i][j+1];
          temp3 = Lele->uMatrix.daa[i][j+2];
          Lele->uMatrix.daa[i][j]   = temp1*cx - temp2*cx*cy/cxz - temp3*cz/cxz;
          Lele->uMatrix.daa[i][j+1] = temp1*cy + temp2*cxz;
          Lele->uMatrix.daa[i][j+2] = temp1*cz - temp2*cy*cz/cxz + temp3*cx/cxz;
       }
    }
    return ( Lele );
}