elmt_shell_4n.c

有限元程序

	     break;

        case MASS_MATRIX:  /* form mass matrix */

#ifdef DEBUG
       printf("*** In elmt_shell_4nodes_implicit() : start case MASS\n");
       printf("                : Density = %14.4e\n", density.value);
       printf("                : shell_thickness = %8.2f\n", h);
       printf("                : Jac             = %8.2f\n", Jac);
#endif
   /*====================================================*/
   /*  CALCULATING MASS MATRIX IN CO_COORDINATE SYSTEM   */ 
   /*====================================================*/

      Shell_4Node_Mass(p, p->stiff, co_coord, density.value, h);

#ifdef DEBUG
       printf("*** In elmt_shell() : end case MASS\n");
#endif
             break;
        default:
             break;
    }

    MatrixFreeIndirectDouble(co_coord, p->no_dimen);
    MatrixFreeIndirectDouble(Direction_Matrix,6);

#ifdef DEBUG
       printf("*** leaving elmt_shell() \n");
#endif

    return(p);
}

/* Print SHELL_4N Element Properties */
#ifdef __STDC__
void print_property_shell_4n(EFRAME *frp, int i)
#else
void print_property_shell_4n(frp, i)
EFRAME    *frp;
int          i;                 /* elmt_attr_no */
#endif
{
int     UNITS_SWITCH;
ELEMENT_ATTR    *eap;

#ifdef DEBUG
       printf("*** Enter print_property_shell_4n()\n");
#endif

     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]);
        } 
     }

     switch(UNITS_SWITCH) {
       case ON:
        UnitsSimplify( eap->work_material[0].dimen );
        UnitsSimplify( eap->work_material[2].dimen );
        UnitsSimplify( eap->work_material[5].dimen );
        UnitsSimplify( eap->work_section[2].dimen );
        UnitsSimplify( eap->work_section[10].dimen );
        if( eap->work_material[0].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Young's Modulus =  E = %16.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[4].value != 0.0 ) {
           printf("             ");
           printf("         : Poisson's ratio = nu = %16.3e   \n", eap->work_material[4].value);
        }
        if( eap->work_material[2].dimen->units_name != NULL ) {
           printf("             ");
           printf("         : Yielding Stress = fy = %16.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[5].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Density         = %16.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_section[2].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Inertia Izz     = %16.3e %s\n",
                           eap->work_section[2].value/eap->work_section[2].dimen->scale_factor,
                           eap->work_section[2].dimen->units_name);
	}
	if( eap->work_section[10].dimen->units_name != NULL ) {
          printf("             ");
          printf("         : Area            = %16.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:
         printf("             ");
         printf("         : Young's Modulus =  E = %16.3e\n",
                            eap->work_material[0].value);
         printf("             ");
         printf("         : Yielding Stress = fy = %16.3e\n",
                            eap->work_material[2].value);
         printf("             ");
         printf("         : Poisson's ratio = nu = %16.3e   \n", eap->work_material[4].value);
         printf("             ");
         printf("         : Density         = %16.3e\n",
                            eap->work_material[5].value);
         printf("             ");
         printf("         : Inertia Izz     = %16.3e\n",
                            eap->work_section[2].value);
         printf("             ");
         printf("         : Area            = %16.3e\n",
                            eap->work_section[10].value);
        break;
        default:
        break;
     }
#ifdef DEBUG
       printf("*** Leave print_property_shell_4n()\n");
#endif
}

/* ============================================*/
/* Equivalent Load due to distributed pressure */
/* for shell Belytschko element                */
/* in local coordinate                         */
/* ============================================*/

ARRAY *sld04(p,task)
ARRAY *p;
int task;
{
ELEMENT_LOADS               *elsptr;
ELOAD_LIB                      *elp;

double                  *nodal_load;
double      Jac, x31, y31, x42, y42; 
double px,py,pz, mx,my,mz, bx,by,bz;

int  i,j,k, ii, jj, kk, n, n1,n2,n3;

#ifdef DEBUG
      printf("**** enter sld04(): \n");
#endif

    /* Initialize total load */

    nodal_load = dVectorAlloc(p->size_of_stiff);

    for(i =1; i <= p->size_of_stiff ;i++)
        nodal_load[i-1] = 0.0;

#ifdef DEBUG
      printf(" **** In sld04(): begin switch()  \n");
#endif

    switch(task){
       case PRESSLD:
#ifdef DEBUG
      printf(" **** In sld04(): case PRESSLD\n");
#endif

       /* ==========================================*/
       /* For distrbuted surface loading:           */
       /* PRESSURE OR TRACTION                      */
       /* nodal equivalent force is calculated by : */
       /* force = integral N^T*traction/pressure dA */
       /* ==========================================*/
       
       /*================================================*/
       /* use one-point integration rule in surface area */
       /* Ni = 1/4.0 [N] = 1/4 [I, I, I, I]              */
       /* [I] = 5x5 unit_matrix                          */
       /* Let traction = [p]5x1                          */
       /* Integ [N]^t*traction dA  = (1/4)*[p]*Jac*4.0   */
       /*                                  [p]           */
       /*                                  [p]           */
       /*                                  [p]           */
       /*                                  [p]           */
       /*================================================*/
         elsptr  =  p->elmt_load_ptr;
         for (j=1; j<= elsptr->no_loads_faces;j++) {
              elp = &elsptr->elib_ptr[j-1];
       
              x31 =  p->coord[0][2].value - p->coord[0][0].value;
              y31 =  p->coord[1][2].value - p->coord[1][0].value;

              x42 =  p->coord[0][3].value - p->coord[0][1].value;
              y42 =  p->coord[1][3].value - p->coord[1][1].value;

              Jac = 0.5*(x31*y42-x42*y31);
 
              for(i = 1; i <= p->size_of_stiff; i++) {
                  n = i/p->dof_per_node;
                  k = i - n*p->dof_per_node;
                  nodal_load[i-1] = 0.25*elp->traction[k-1].value*4.0*Jac; 
              }
         }
       break;
    default:
       break;
    }

    for(i = 1; i <= p->size_of_stiff; i++) {
       p->equiv_nodal_load->uMatrix.daa[i-1][0] = nodal_load[i-1];
    }
    
    free(nodal_load);

#ifdef DEBUG
      printf("**** leaving sld04(): \n");
#endif
    return(p);
}


/* ============================*/
/* Material matrix             */
/* ============================*/

double **MATER_MAT_SHELL(m1, E, nu)
double  **m1;
double E, nu;
{
double           temp;
double       kk = 1.2;       
    
#ifdef DEBUG
     printf(" enter MATER_MAT_SHELL() \n");
     printf(" E = %lf , nu = %lf \n", E, nu);
#endif
    /* Stress : stress_x, stress_y,  stress_xy,  stress_yz,  stress_zx  */
    /* Strain : strain_x, strain_y, 2strain_xy, 2strain_yz, 2strain_zx  */

    
    temp =  E/(1-nu*nu);

    m1[0][0]= m1[1][1] = temp;
    m1[0][1]= m1[1][0] = nu*temp;
    m1[2][2]= E/2.0/(1.0+nu);
    m1[3][3] = m1[4][4] = E/2.0/(1.0+nu)/kk;

    m1[0][2] = m1[2][0] = m1[1][2] = m1[2][1] = 0.0;
    m1[3][0] = m1[3][1] = m1[3][2] = m1[3][4] = 0.0;
    m1[4][0] = m1[4][1] = m1[4][2] = m1[4][3] = 0.0;


#ifdef DEBUG
     dMatrixPrint("m1", m1, 5, 5);
     printf("\n leaving MATER_MAT_SHELL() \n");
#endif
  return (m1);
}

void MATER_SHELL_UPDATE(p, co_coord, E, nu, integ_pt)
ARRAY                              *p;
double                     **co_coord;
double                          E, nu;
int                          integ_pt; /* integration point */
{
double             **m1, **mater_temp;
double **stress_dev, **Norm_Transpose;
double         A_2, mean_stress, R, A;
double                           temp;
int            i, j, k, ii, jj, kk, n;
int                      dof, size, G;

#ifdef DEBUG
    printf(" enter MATER_SHELL_UPDATE() \n");
#endif

  /* Given i -state displ, stress, [C] and strain */
  /* calculate i+1 state displ, stress, these i+1 */
  /* state stress can be used for i+1 state [Cep] */
  /* matrix                                       */

#ifdef DEBUG
    printf(" In MATER_SHELL_UPDATE(): calculating MATER_MAT_SHELL() \n");
#endif
      dof             = p->dof_per_node;
      size            = p->size_of_stiff;
      ii              = integ_pt;
      G               = E/2.0/(1.0+nu);


#ifdef DEBUG
    printf(" In MATER_SHELL_UPDATE(): update p->material = %s \n", p->material_name);
    printf(" In MATER_SHELL_UPDATE(): integ_pts          = %d \n", integ_pt);
    printf(" In MATER_SHELL_UPDATE(): dof                = %d \n", dof);
    printf(" In MATER_SHELL_UPDATE(): deformation state  = %d \n", p->elmt_state);
#endif

  /* Elastic material */

  if(p->material_name != NULL && !strcmp(p->material_name, "ELASTIC")) {

     p->mater_matrix->uMatrix.daa 
     = MATER_MAT_SHELL(p->mater_matrix->uMatrix.daa,E, nu);
   }

  /* Elastic Perfectly Plastic Materials */


  if(p->material_name != NULL &&
     !strcmp(p->material_name,"ELASTIC_PERFECTLY_PLASTIC")) {

     m1              = MatrixAllocIndirectDouble(p->dof_per_node, p->dof_per_node);
     stress_dev      = MatrixAllocIndirectDouble(dof, 1);
     Norm_Transpose  = MatrixAllocIndirectDouble(1,dof);
     mater_temp      = MatrixAllocIndirectDouble(dof, dof);

     /* Elastic Deformation */

     switch(p->elmt_state) {

       /* Elastic deformation */
       case 0: 
         p->mater_matrix->uMatrix.daa 
         = MATER_MAT_SHELL(p->mater_matrix->uMatrix.daa,E, nu);
       break;
       case 1:
       /* plastic deformation */

       /* Step 1: Calculate the spherical stress and */
       /*         deviatric stress of trial stress   */
       /*         also the radiaus square A_2        */

            mean_stress = (p->stress->uMatrix.daa[0][ii-1] - p->LC_ptr->back_stress[0][ii-1]
                          +p->stress->uMatrix.daa[1][ii-1] - p->LC_ptr->back_stress[1][ii-1])/3.0;
            A_2 = 0.0;
            for(i = 1; i <= dof; i++) {
               if(i <= 2)           
                 stress_dev[i-1][0] = p->stress->uMatrix.daa[i-1][ii-1] - mean_stress - p->LC_ptr->back_stress[i-1][ii-1];
               else