elmt_fiber3d.c

有限元程序

        if( eap->work_material[4].value != 0.0 ) {
           printf("             ");
           printf("         : Poisson's Ratio         = nu = %11.3e   \n", eap->work_material[4].value);
        }
        if( eap->work_material[0].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Young's Modulus         =  E = %11.3e %s\n",
                           eap->work_material[0].value/eap->work_material[0].dimen->scale_factor,
                           eap->work_material[0].dimen->units_name);
        }
        if( eap->work_material[3].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Young's Tangent Modulus = Et = %11.3e %s\n",
                           eap->work_material[3].value/eap->work_material[3].dimen->scale_factor,
                           eap->work_material[3].dimen->units_name);
        }
        if( eap->work_material[2].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Yielding Stress         = fy = %11.3e %s\n",
                           eap->work_material[2].value/eap->work_material[2].dimen->scale_factor,
                           eap->work_material[2].dimen->units_name);
        }
        if( eap->work_material[1].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Shear Modulus           =  G = %11.3e %s\n",
                           eap->work_material[1].value/eap->work_material[1].dimen->scale_factor,
                           eap->work_material[1].dimen->units_name);
        }
        if( eap->work_material[10].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Shear Tangent Modulus   = Gt = %11.3e %s\n",
                           eap->work_material[10].value/eap->work_material[10].dimen->scale_factor,
                           eap->work_material[10].dimen->units_name);
        }
        if( eap->work_material[11].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Shear Yielding Stress   = fv = %11.3e %s\n",
                           eap->work_material[11].value/eap->work_material[11].dimen->scale_factor,
                           eap->work_material[11].dimen->units_name);
        }
        if( eap->work_section[16].value != 0.0 ) {
           printf("             ");
           printf("         : Shear Correction Factor = ks = %11.3e   \n", eap->work_section[16].value);
        }
	if( eap->work_section[10].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Section Area            =      %11.3e %s\n",
                           eap->work_section[10].value/eap->work_section[10].dimen->scale_factor,
                           eap->work_section[10].dimen->units_name);
	}
       break;
       case OFF:
	if( eap->work_material[5].value != 0.0 ) {
          printf("             ");
          printf("         : Density                 =      %11.3e\n", eap->work_material[5].value);
	}
        if( eap->work_material[4].value != 0.0 ) {
           printf("             ");
           printf("         : Poisson's Ratio         = nu = %11.3e\n", eap->work_material[4].value);
        }
        if( eap->work_material[0].value != 0.0 ) {
           printf("             ");
           printf("         : Young's Modulus         =  E = %11.3e\n", eap->work_material[0].value);
        }
        if( eap->work_material[3].value != 0.0 ) {
           printf("             ");
           printf("         : Young's Tangent Modulus = Et = %11.3e\n", eap->work_material[3].value);
        }
        if( eap->work_material[2].value != 0.0 ) {
           printf("             ");
           printf("         : Yielding Stress         = fy = %11.3e\n", eap->work_material[2].value);
        }
        if( eap->work_material[1].value != 0.0 ) {
           printf("             ");
           printf("         : Shear Modulus           =  G = %11.3e\n", eap->work_material[1].value);
        }
        if( eap->work_material[10].value != 0.0 ) {
           printf("             ");
           printf("         : Shear Tangent Modulus   = Gt = %11.3e\n", eap->work_material[10].value);
        }
        if( eap->work_material[11].value != 0.0 ) {
           printf("             ");
           printf("         : Shear Yielding Stress   = fv = %11.3e\n", eap->work_material[11].value);
        }
        if( eap->work_section[16].value != 0.0 ) {
           printf("             ");
           printf("         : Shear Correction Factor = ks = %11.3e   \n", eap->work_section[16].value);
        }
	if( eap->work_section[10].value != 0.0 ) {
          printf("             ");
          printf("         : Section Area            =      %11.3e\n", eap->work_section[10].value);
	}
        break;
        default:
        break;
     }

     /* Print fiber attribution */

     printf("             ");
     printf("Fiber Attribution\n");
     printf("             ");
     printf("         : No. of Fibers = %6i\n", eap->work_fiber->no_fiber );
}

/*
 *  =====================================
 *  Fiber_Elmt_State_Det_3d()
 *
 *  Input :
 *  Output : void 
 *  =====================================
 */

#ifdef  __STDC__
void Fiber_Elmt_State_Det_3d( ARRAY *p, HISTORY_DATA *hp, int *flag )
#else
void Fiber_Elmt_State_Det_3d( p , hp, flag)
ARRAY *p;
HISTORY_DATA *hp;
int *flag;
#endif
{
int         no_integ_pt, no_section, no_fiber, no_shear, elmt_no;
int         total_fiber;
QUANTITY    length;

int         ii, jj, kk, sec, ifib;
int         UNITS_SWITCH, UnitsType;
double      cs, sn, tn;
MATRIX      *temp_m1, *temp_m2;
DIMENSIONS  *dp_length, *dp_force, *dp_stress, *dimen;

double      *abscissas, *weights;
double      xi;
double      scale;
MATRIX      *Q, *q;
MATRIX      *dpe, *dQ, *dq, *s;
MATRIX      *Dx, *dx, *ex, *rx;
MATRIX      *dDx, *ddx, *dex;
MATRIX      *DRx, *DUx;
MATRIX      *F, *K;
MATRIX      *L, *Ltrans, *R, *Rtrans;
MATRIX      *kx, *fx;
MATRIX      *lx, *bx, *bxtrans;
MATRIX      *stress, *tangent;
FIBER_ATTR  *fiber;
SECTION_DATA *saved_data;

    UNITS_SWITCH = CheckUnits();
    UnitsType    = CheckUnitsType();

    /* Assign Element Property */

    fiber       = p->fiber_ptr->fiber;
    no_fiber    = p->fiber_ptr->no_fiber;
    no_shear    = p->fiber_ptr->no_shear;
    total_fiber = no_fiber + (p->no_dimen-1)*no_shear;
    no_integ_pt = p->integ_ptr->integ_pts;
    no_section  = no_integ_pt + 2;
    elmt_no     = p->elmt_no;
    cs = p->coord[0][1].value - p->coord[0][0].value;
    sn = p->coord[1][1].value - p->coord[1][0].value;
    tn = p->coord[2][1].value - p->coord[2][0].value;
    length.value = sqrt( cs*cs + sn*sn + tn*tn );
    if( UNITS_SWITCH == ON ) {
       if( UnitsType == SI )
          dp_length = DefaultUnits("m");
       else if( UnitsType == US )
          dp_length = DefaultUnits("in");
       length.dimen = (DIMENSIONS *)MyCalloc(1, sizeof(DIMENSIONS));
       UnitsCopy( length.dimen, dp_length );
       free((char *) dp_length->units_name);
       free((char *) dp_length);

       SetUnitsOff();
    }

    /* Gauss-Lobatto Integration */

    abscissas = (double *)MyCalloc( no_section, sizeof(double) );
    weights   = (double *)MyCalloc( no_section, sizeof(double) );
    Gauss_Lobatto( abscissas, weights, no_integ_pt );

    /* Put Displacement Matrix To One Column Form */
    /* p->displ=[dpx1,dpx2; dpy1,dpy2; dpz1,dpz2; drx1,drx2; dry1,dry2; drz1,drz2]6x2 */
    /*  =>  dpe=[dpx1;dpy1;dpz1;drx1;dry1;drz1; dpx2;dpy2;dpz2;drx2;dry2;drz2]12x1   */

    dpe = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 12, 1 );
    for( ii=0 ; ii < p->dof_per_node ; ++ii ) {
       for( jj=0 ; jj < p->nodes_per_elmt ; ++jj ) {
          kk = ii + jj*p->dof_per_node;
          dpe->uMatrix.daa[kk][0] = p->displ->uMatrix.daa[ii][jj];
       }
    }

    /* Use the same address in array, therefore  */
    /* frame->element[elmt_no-1]->rp->Q_saved, q_saved will update automatically. */
    Q = p->Q_saved;
    q = p->q_saved;
    ex = hp->strain;   /* ex[no_section][total_fiber] */
    stress = hp->stress;
    tangent = hp->tangent;
#ifdef DEBUG
printf("\n Q_saved\n");
MatrixPrintVar( Q, (MATRIX *)NULL );
printf("\n q_saved\n");
MatrixPrintVar( q, (MATRIX *)NULL );
printf("\n stress\n");
MatrixPrintVar( stress, (MATRIX *)NULL );
printf("\n tangent\n");
MatrixPrintVar( tangent, (MATRIX *)NULL );
#endif

    /* Calculate The Initial Tangent Stiffness, And Each Section Related Matrix */

    saved_data = (SECTION_DATA *)MyCalloc( no_section, sizeof(SECTION_DATA) );
    kx = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 5, 5 );
    bx = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 5, 5 );
    lx = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, total_fiber, 5 );

    for( sec=0 ; sec < no_section; ++sec ) {
       saved_data[sec].xi = length.value/2.0*(abscissas[sec]+1);
       saved_data[sec].wi = weights[sec]*length.value/2.0;

       Force_Interpolation_Matrix_3d( bx, saved_data[sec].xi, length );
       Section_Tangent_Stiffness_3d( kx, fiber, no_fiber, no_shear, sec, elmt_no );

       saved_data[sec].fx = MatrixInverse( kx );
       saved_data[sec].Dx = MatrixMult( bx, Q );
       dx = MatrixMult( saved_data[sec].fx, saved_data[sec].Dx );
       saved_data[sec].rx = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 5, 1 );
       MatrixFree( dx );

       bxtrans = MatrixTranspose( bx );
       temp_m1 = MatrixMult( bxtrans, saved_data[sec].fx );
       temp_m2 = MatrixMult( temp_m1, bx );
       MatrixFree( temp_m1 );
       MatrixFree( bxtrans );
       if( sec == 0 )
          F = MatrixScale( temp_m2, saved_data[sec].wi );
       else {
          temp_m1 = MatrixScale( temp_m2, saved_data[sec].wi );
          MatrixAddReplace( F, temp_m1 );
          MatrixFree( temp_m1 );
       } 
       MatrixFree( temp_m2 );
    }
    K = MatrixInverse( F );

    /* Element Deformation Increments */

    L = Element_Transformation_3d( p->coord, length );
    dq = MatrixMult( L, dpe );
    MatrixAddReplace( q, dq );
#ifdef DEBUG
printf("\n dpe\n");
MatrixPrintVar( dpe, (MATRIX *)NULL );
printf("\n dq\n");
MatrixPrintVar( dq, (MATRIX *)NULL );
#endif

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

    s   = MatrixAllocIndirect(   "s",  DOUBLE_ARRAY, 5, 1 );
    DRx = MatrixAllocIndirect( "DRx",  DOUBLE_ARRAY, 5, 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 DEBUG
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 DEBUG
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_3d( bx, xi, length );
          Linear_Geometric_Matrix_3d( 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 DEBUG
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] );