elmt_fiber2d.c

有限元程序

                FatalError("FIBER_2D element only support LUMPED mass", (char *)NULL);
                break;
          }

          break;

       default:
          printf("*** ERROR in elmt_fiber2d() : isw = %d not defined\n", isw );
          break;
    }

    return(p);
}


/* ========================================================== */
/* print_property_fiber_2d() : Print Fiber Element Properties */
/*                                                            */
/* Input:                                                     */
/*                                                            */
/*   EFRAME *frp :                                            */
/*   int i       :                                            */
/*                                                            */
/* Output: void.                                              */
/* ========================================================== */

#ifdef __STDC__
void print_property_fiber_2d(EFRAME *frp, int i)
#else
void print_property_fiber_2d(frp, i)
EFRAME    *frp;
int          i;                 /* elmt_attr_no */
#endif
{
int     UNITS_SWITCH;
ELEMENT_ATTR    *eap;
int             ifib;

     UNITS_SWITCH = CheckUnits();
     eap = &frp->eattr[i-1];

     if( PRINT_MAP_DOF == ON ) {
        if(frp->no_dof == 3 || frp->no_dof == 2) { 
           printf("             ");
           printf("         : gdof [0] = %4d : gdof[1] = %4d : gdof[2] = %4d\n",
                           eap->map_ldof_to_gdof[0],
                           eap->map_ldof_to_gdof[1],
                           eap->map_ldof_to_gdof[2]);
        }

        if(frp->no_dof == 6) { /* 3d analysis */
           printf("             ");
           printf("         : dof-mapping : gdof[0] = %4d : gdof[1] = %4d : gdof[2] = %4d\n",
                           eap->map_ldof_to_gdof[0],
                           eap->map_ldof_to_gdof[1],
                           eap->map_ldof_to_gdof[2]);
           printf("             ");
           printf("                         gdof[3] = %4d : gdof[4] = %4d : gdof[5] = %4d\n",
                           eap->map_ldof_to_gdof[3],
                           eap->map_ldof_to_gdof[4],
                           eap->map_ldof_to_gdof[5]);
        } 
     }

     printf("             ");
     printf(" General Section and Material Properties\n");
     switch(UNITS_SWITCH) {
       case ON:
        UnitsSimplify( eap->work_material[0].dimen );
        UnitsSimplify( eap->work_material[1].dimen );
        UnitsSimplify( eap->work_material[2].dimen );
        UnitsSimplify( eap->work_material[3].dimen );
        UnitsSimplify( eap->work_material[5].dimen );
        UnitsSimplify( eap->work_material[10].dimen );
        UnitsSimplify( eap->work_material[11].dimen );
        UnitsSimplify( eap->work_section[10].dimen );
	if( eap->work_material[5].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Density                 =      %11.3e %s\n",
                 eap->work_material[5].value/eap->work_material[5].dimen->scale_factor,
                 eap->work_material[5].dimen->units_name);
	}
        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 Attributes1 */

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


/*
 *  =================================================================
 *  State Determination for 2D fiber element
 *
 *  Input  : Array        *p  -- pointer to working array.
 *         : HISTORY_DATA *pp -- pointer to history info.
 *         : int *flag        -- flag.
 *  Output : void
 *  =================================================================
 */

/* #define DEBUG2 */

#ifdef  __STDC__
void Fiber_Elmt_State_Det_2d( ARRAY *p, HISTORY_DATA *hp, int *flag )
#else
void Fiber_Elmt_State_Det_2d( 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, xx, yy;
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;

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

    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;
    length.value = sqrt( cs*cs + sn*sn );
    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]3x2 */
    /*  =>  dpe=[dpx1;dpy1;dpz1;dpx2;dpy2;dpz2]6x1   */

    dpe = MatrixAllocIndirect( (char *)NULL, DOUBLE_ARRAY, 6, 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 DEBUG2
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 related section matrix */

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

    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_2d( bx, saved_data[sec].xi, length );
       Section_Tangent_Stiffness_2d( 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, 3, 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_2d( p->coord, length );
    dq = MatrixMult( L, dpe );
    MatrixAddReplace( q, dq );

#ifdef DEBUG2
printf("\n dpe\n");
MatrixPrintVar( dpe, (MATRIX *)NULL );
printf("\n dq\n");
MatrixPrintVar( dq, (MATRIX *)NULL );