elmt_psps.c

有限元程序

               for(k = 1; k <= nodes_per_elmt; k++)
                   xs[i-1][j-1] = xs[i-1][j-1] + coord[i-1][k-1].value*shp_temp[j-1][k-1];
           }

           *jacobian = xs[0][0] * xs[1][1] - xs[0][1] *xs[1][0];

           if(flg == TRUE)
              return;

           /* Compute Jacobian inverse matrix */

           sx[0][0] = xs[1][1]/ *jacobian;
           sx[1][1] = xs[0][0]/ *jacobian;
           sx[0][1] = - xs[0][1]/ *jacobian;
           sx[1][0] = - xs[1][0]/ *jacobian;

           /* Save dNi/dx, dNi/dy into shp[2][node] */
  
           for(i = 1; i <= nodes_per_elmt; i++){
               shp[0][i-1] = shp_temp[0][i-1]*sx[0][0] + shp_temp[1][i-1]*sx[0][1];
               shp[1][i-1] = shp_temp[0][i-1]*sx[1][0] + shp_temp[1][i-1]*sx[1][1];
           }
           break;
      default:
           break;
   }

   MatrixFreeIndirectDouble(shp_temp, 2);
   return;
}


/*
 *  ============================================================================= 
 *  gauss() : Gauss Integration 
 *
 *  Input :
 *  Output :
 *  ============================================================================= 
 */

#ifdef __STDC__
gauss( double *sg, double *ag, int lt)
#else
gauss(sg, ag, lt)
double *sg, *ag;
int lt;
#endif
{
double t;
int    i;

    switch(lt) {
        case 1:
             sg[1] = 0.0;
             ag[1] = 2.0;
             break;
        case 2:
             sg[1] = -1/sqrt(3.);
             sg[2] =  -sg[1];
             ag[1] = 1.0; 
             ag[2] = 1.0; 
             break;
        case 3:
             sg[1] = -sqrt(0.6);
             sg[2] = 0.0;
             sg[3] = -sg[1];

             ag[1] = 5.0/9.0; 
             ag[2] = 8.0/9.0; 
             ag[3] = ag[1]; 
             break;
        case 4:
             t = sqrt(4.8);
             sg[1] =  sqrt((3+t)/7);
             sg[2] =  sqrt((3-t)/7);
             sg[3] = -sg[2];
             sg[4] = -sg[1];

             t  = 1/3/t;
             ag[1] = 0.5-t; 
             ag[2] = 0.5+t; 
             ag[3] = ag[2]; 
             ag[4] = ag[1]; 
             break;
        default:
             break;
    }

    return (1);
}

/*
 *  ============================================================================= 
 *  pgauss() : 
 * 
 *  Input  : int no_integ_pts --  no of gaussian integ pts in one-direction.
 *  Output : lint             --  no_integ_pts*no_integ_pts  
 *         :   sg             -- coordinates in xi direction 
 *         :   tg             -- coordinates in eta direction
 *         :   wg             -- weighting coefficients
 *  ============================================================================= 
 */

int pgauss(l, lint, r, z, w)
int l, *lint;
double r[], z[], w[];
{
int i, j, k;
double g, h;
double g4[5], h4[5], lr[10], lz[10], lw[10];

lr[1] = lr[4] =lr[8] = -1;
lr[2] = lr[3] =lr[6] =  1;
lr[5] = lr[7] =lr[9] =  0;

lz[1] = lz[2] =lz[5] = -1;
lz[3] = lz[4] =lz[7] =  1;
lz[6] = lz[8] =lz[9] =  0;

lw[3] = lw[4] = lw[1] = lw[2] = 25;
lw[5] = lw[6] = lw[7] = lw[8] = 40;
lw[9] = 64;

    if(l < 0) {
       *lint = 3;
       g = sqrt((double) 1.0/3.0);
       h = sqrt((double) 2.0/3.0);
       r[1] = -h; r[2] =  h; r[3] = 0;
       z[1] = -g; z[2] = -g; z[3] = g;
       w[1] =  1; w[2] =  1; w[3] = 2;

       return (1);
    }

    *lint = l*l;
    switch (l) {
	case 1: /*  1 x 1 integration */
	     r[0] = 0.0;
	     z[0] = 0.0;
	     w[0] = 4.0;
	     break;
        case 2: /*  2 x 2 integration */
             g = 1.0/sqrt((double) 3.0);
             for(i = 1; i<= 4; i++) {
                 r[i-1] = g * lr[i];
                 z[i-1] = g * lz[i];
                 w[i-1] = 1.0;
             }
             break;
        case 3: /* 3 x 3 integration */
             g = sqrt((double) 0.6);
             h = 1.0/81.0;
             for(i = 1; i<= 9; i++) {
                 r[i-1] = g * lr[i];
                 z[i-1] = g * lz[i];
                 w[i-1] = h * lw[i];
             }
             break;
        case 4: /* 4 x 4 integration */
             g = sqrt((double) 4.8);
             h = sqrt((double) 30.0)/36;
             g4[1] = sqrt((double) (3+g)/7.0);
             g4[4] = -g4[1];
             g4[2] = sqrt((double) (3-g)/7.0);
             g4[3] = -g4[2];
             h4[1] = 0.5 - h;
             h4[2] = 0.5 + h;
             h4[3] = 0.5 + h;
             h4[4] = 0.5 - h;
             i = 0;
             for(j = 1; j<= 4; j++) {
                 for(k = 1; k<= 4; k++) {
                     i = i +1;
                     r[i-1] = g4[k];
                     z[i-1] = g4[j];
                     w[i-1] = h4[j]* h4[k];
                 }
             }
             break;
    }

    return(1);
}

/*
 *  ===============================================================
 *  MaterialMatrixPlane() : Compute (3x3) constituitive matrix.
 *
 *  Note : Plane Stress : dFlag = 1
 *         Plane Strain : dFlag = 2
 *
 *  Input  : double     E -- Moduus of elasticity.
 *         : double    nu -- Poisson's ratio
 *         : double dFlag -- Flag for Plane Strss/Plane Strain Cases
 *  Output : double **m1  -- Pointer to constituitive matrix.
 *  ===============================================================
 */

double **MaterialMatrixPlane( double E, double nu , double dFlag ) {
double **m1;
double temp;

    m1 = MatrixAllocIndirectDouble(3, 3);
    
    temp =  E*(1 + ( 1-dFlag)*nu)/(1.0+nu)/(1.0 - dFlag*nu);

    m1[0][0] = m1[1][1] = temp;
    m1[0][1] = m1[1][0] = (nu/(1+(1-dFlag)*nu))*temp; 
    m1[2][2] = E/2.0/(1.0+nu); 

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

    return(m1);
}


/*
 *  ==========================================================================
 *  print_property_psps() : print PLANE_STRAIN/PLANE_STRESS element properties
 *  ==========================================================================
 */

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

     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:
        if( eap->work_material[0].value != 0.0 ) {
           printf("             ");
           printf("         : Young's Modulus =  E = %16.3e\n",
                            eap->work_material[0].value);
        }
        if( eap->work_material[2].value != 0.0 ) {
           printf("             ");
           printf("         : Yielding Stress = fy = %16.3e\n",
                            eap->work_material[2].value);
        }
        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[0].value != 0.0 ) {
           printf("             ");
           printf("         : Density         = %16.3e\n",
                            eap->work_material[5].value);
        }
        if( eap->work_section[2].value != 0.0 ) {
           printf("             ");
           printf("         : Inertia Izz     = %16.3e\n",
                            eap->work_section[2].value);
        }
        if( eap->work_section[10].value != 0.0 ) {
           printf("             ");
           printf("         : Area            = %16.3e\n",
                            eap->work_section[10].value);
        }
        break;
        default:
        break;
     }
}