elmt_shell_8n.c

有限元程序

    v3_ptr[0][0] = s_ptr[0][0]/s_norm;
    v3_ptr[1][0] = s_ptr[1][0]/s_norm;
    v3_ptr[2][0] = s_ptr[2][0]/s_norm;

#ifdef DEBUG
    dMatrixPrint("v3_ptr", v3_ptr, 3,1);
#endif

 /*     v1_ptr = ey_ptr X v3_ptr      */

    v1_ptr = dVmatrixCrossProduct(v1_ptr, ey_ptr, 3, 1, v3_ptr, 3,1);

    s_norm  = (double) dVmatrixL2Norm(v1_ptr, 3, 1);

    /* check if ey_ptr is parallel to v3_ptr */

    if(abs(s_norm) <1E-7) { /* parallel */
       v1_ptr = dVmatrixCrossProduct(v1_ptr, ez_ptr, 3, 1, v3_ptr, 3,1);
    }

    v2_ptr = dVmatrixCrossProduct(v2_ptr, v3_ptr, 3, 1, v1_ptr, 3,1);

    /* assume that each node has the same directions in */
    /* one elements                                     */

    for (i = 1; i <= 3; i++) {
         for (k = 1; k <= p->nodes_per_elmt; k++) {
             e1_ptr[i-1][k-1] = v1_ptr[i-1][0];
             e2_ptr[i-1][k-1] = v2_ptr[i-1][0];
             e3_ptr[i-1][k-1] = v3_ptr[i-1][0];
        }
    }

    MatrixFreeIndirectDouble(s_ptr, 3);

    MatrixFreeIndirectDouble(v1_ptr, 3);
    MatrixFreeIndirectDouble(v2_ptr, 3);
    MatrixFreeIndirectDouble(v3_ptr, 3);

    MatrixFreeIndirectDouble(ez_ptr, 3);
    MatrixFreeIndirectDouble(ey_ptr, 3);

    MatrixFreeIndirectDouble(r75_ptr, 3);
    MatrixFreeIndirectDouble(r86_ptr, 3);

#ifdef DEBUG
        printf(" Leaving Lamina_Sys_8node()\n");
#endif

}

/* ======================*/
/* Shell Stiff Matrix    */
/* ======================*/

#ifdef __STDC__
void Shell_Stiff_Plane_8node(double **K, ARRAY *p, double **co_coord, double z_coord, int z_integ_pt, double EE, double nu)
#else
void Shell_Stiff_Plane_8node(K, p, co_coord, z_coord, z_integ_pt, EE, nu)
double         **K;
ARRAY           *p;
double  **co_coord;
double     z_coord;      /* gussain pt in z-direction */  
int     z_integ_pt;      /* integration point No in z-direction */
double      EE, nu;
#endif
{
int                i, j, k,n, ii, jj, kk;
static int     surface_pts, no_integ_pts;
int                            size, dof;
double                    *x_integ_coord;
double                    *y_integ_coord;
double                 **J_inverse= NULL;
static double       weight[16], jacobian;
static double               **shp = NULL;
static double      **stiff_matrix = NULL;
static double                **TT = NULL;
static double          **B_matrix = NULL;
static double         **B1_matrix = NULL; /* local B_matrix           */
static double         **B1_Trans  = NULL; /* local B_matrix Transpose */
static double       **temp_matrix = NULL;

static double    diff, sum, temp1, temp2;

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

    shp            = MatrixAllocIndirectDouble(3,8);
    J_inverse      = MatrixAllocIndirectDouble(3,3);
    TT             = MatrixAllocIndirectDouble(6,6);
    x_integ_coord  = dVectorAlloc(16);
    y_integ_coord  = dVectorAlloc(16);

    stiff_matrix   = MatrixAllocIndirectDouble(p->size_of_stiff,p->size_of_stiff);
    B_matrix       = MatrixAllocIndirectDouble(p->dof_per_node+ 1,p->size_of_stiff);
    B1_matrix      = MatrixAllocIndirectDouble(p->dof_per_node, p->size_of_stiff);
    B1_Trans       = MatrixAllocIndirectDouble(p->size_of_stiff,p->dof_per_node);
    temp_matrix    = MatrixAllocIndirectDouble(p->dof_per_node,p->size_of_stiff);

    surface_pts    = p->integ_ptr->surface_pts;
    no_integ_pts   = (int) 0;
    size           = p->size_of_stiff;
    dof            = p->dof_per_node;

    if(surface_pts*surface_pts != no_integ_pts)
       pgauss(surface_pts, &no_integ_pts, x_integ_coord, y_integ_coord, weight);

    /* Material Matrix for elastic materials  */

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

   /* ==============================*/
   /* initialize stiffness matrix   */
   /* ==============================*/
    for ( i = 1; i <= size; i++) {
         for(j = 1; j <= size; j++) {
             K[i-1][j-1] = 0.0;
         }
    }

    for( ii = 1; ii <= no_integ_pts; ii++) { /* loop over the surface */

        elmt_shell_shape_8node(p, co_coord, x_integ_coord[ii-1],
                        y_integ_coord[ii-1], z_coord, shp,&jacobian,J_inverse, TT, STIFF);

        B_matrix = B_MATRIX_8Node(B_matrix, p, shp, z_coord, J_inverse);

       /* Calculate the B_matrix in local coordinate */
       /* eliminate the 3rd row of TT matrix         */

       for(i = 1; i <= 2; i++) {
           for(j = 1; j <= size; j++) {
               B1_matrix[i-1][j-1] = 0.0;
               B1_Trans[j-1][i-1]  = 0.0;
               for(k = 1; k <= dof + 1; k++)
                   B1_matrix[i-1][j-1] += TT[i-1][k-1]*B_matrix[k-1][j-1];
           }
       }

       for(i = 3; i <= dof; i++) { 
           for(j = 1; j <= size; j++) {
               B1_matrix[i-1][j-1] = 0.0;
               B1_Trans[j-1][i-1]  = 0.0;
               for(k = 1; k <= dof + 1; k++)
                   B1_matrix[i-1][j-1] += TT[i][k-1]*B_matrix[k-1][j-1];
           }
       }

#ifdef DEBUG
       printf(" jacobian = %lf \n", jacobian);
       dMatrixPrint(" shape func", shp,3,8);
       dMatrixPrint(" B_matrix", B_matrix, p->dof_per_node + 1, p->size_of_stiff);
#endif
       for( i = 1; i <= dof; i++) {
           for(j = 1; j <= size; j++) {
               B1_Trans[j-1][i-1] = B1_matrix[i-1][j-1];
           }
       }

       /* ============================================================= */
       /* Update the Material Matrix for Elastic_Plastic Materials      */
       /* ============================================================= */

       kk = no_integ_pts*(z_integ_pt-1) +ii;
       MATER_SHELL_UPDATE(p, co_coord, EE, nu, kk);

       /* calculate the element stiffness matrix */

       temp_matrix  = dMatrixMultRep(temp_matrix,
                      p->mater_matrix->uMatrix.daa, dof, dof, B1_matrix, dof, size); 
       stiff_matrix = dMatrixMultRep(stiff_matrix, B1_Trans,
                      size, dof, temp_matrix, dof, size); 

       for ( i = 1; i <= size; i++) {
            for(j = 1; j<= size; j++) {
                K[i-1][j-1] += stiff_matrix[i-1][j-1]*jacobian*weight[ii-1];
            }
       }

    }  /* end of gaussian integration */

#ifdef DEBUG
               /* check the symmetry of the stiffness matrix */
       printf(" INSIDE CHECK \n");
                for(i = 1; i <= size; i++){
                    for(j = 1; j <= size; j++) {
                        diff = ABS(K[i-1][j-1] - K[j-1][i-1]);

                       if(diff > 1E-7 && ABS(diff/K[i-1][j-1]) > 1E-1 &&
                         ABS(K[i-1][j-1]) > 1E-7) {
                           printf("K[%d][%d] = %le \n", i, j, K[i-1][j-1]);
                           printf("K[%d][%d] = %le \n", j, i, K[j-1][i-1]);
                           printf("diff = %le (diff/K[%d][%d]) = %le\n", diff, i, j, (diff/ABS(K[i-1][j-1])));
                           printf("elmtNo = %d Stiffness matrix IS NOT SYMMETRIC \n", p->elmt_no);
                           break;
                       }
                       else {
                           if( i == size && j == size)
                              printf("elmtNo = %d Stiffness matrix IS SYMMETRIC \n", p->elmt_no);
                       }
                    }
                }
#endif



#ifdef DEBUG
       printf(" end of integration over surface \n"); 
#endif

    free((char *) x_integ_coord);
    free((char *) y_integ_coord);
    MatrixFreeIndirectDouble(shp, 3);
    MatrixFreeIndirectDouble(B_matrix, p->dof_per_node+1);
    MatrixFreeIndirectDouble(B1_matrix, p->dof_per_node);
    MatrixFreeIndirectDouble(B1_Trans, p->size_of_stiff);
    MatrixFreeIndirectDouble(temp_matrix, p->dof_per_node);
    MatrixFreeIndirectDouble(stiff_matrix, p->size_of_stiff);

#ifdef DEBUG 
     dMatrixPrint("stiff surface ", K, p->size_of_stiff, p->size_of_stiff);
     printf(" Leaving Shell_Stiff_Plane_8node() \n");
#endif

}


#define MASS    1

/* =================== */
/* Shell Mass Matrix    */
/* =================== */
#ifdef __STDC__
void  Shell_8Node_Mass(ARRAY *p, MATRIX *mass, double **co_coord, double density, double thickness)
#else
void  Shell_8Node_Mass(p, mass, co_coord, density, thickness)
ARRAY               *p;
MATRIX           *mass;
double      **co_coord;
double         density;
double       thickness;
#endif
{
int  i, j, k, n, ii,jj, length, length1, length2;
int    no_integ_pts, x_integ_pts, y_integ_pts;
int                               z_integ_pts;

double                     **J_inverse = NULL;
static double                     **TT = NULL;
static double                 **e1_ptr = NULL;
static double                 **e2_ptr = NULL;
static double                 **e3_ptr = NULL;

double             **shp, **Mass, **Mass_temp;
double                      **N1_Trans = NULL;
double                      **N2_Trans = NULL;
double               **N1 = NULL, **N2 = NULL;

double         *x_integ_coord, *y_integ_coord;
double                  *z_integ_coord = NULL;
double                       *z_weight = NULL;
double                   weight[16], jacobian;
double                                    sum;

double               Aera, elmt_length, width;
double            node_mass, node_Jx, node_Jy;

DIMENSIONS                      *d1, *d2, *d3;
int                   UNITS_SWITCH, UnitsType;

#ifdef DEBUG
       printf("*** Enter Shell_8Node_Mass(): \n");
#endif
     /* [a] mass initialization */

      for(i = 1; i <= p->size_of_stiff; i++) 
          for(j = 1; j <= p->size_of_stiff; j++)
                mass->uMatrix.daa[i-1][j-1] = 0.0;

    /* [b] :  mass matrix      */

#ifdef DEBUG
       printf("*** In Shell_8Node_Mass(): in main swicth() \n");
       printf("*** mtype = %d \n", p->type);
       printf("*** density   = %lf\n", density);
       printf("*** thickness = %lf\n", thickness);
#endif


    switch(p->type) {
	case LUMPED:
              
             for (i = 1; i <= p->nodes_per_elmt; i++) { 
                  elmt_length = ABS(p->coord[0][i-1].value - p->coord[0][i].value);
                  if(elmt_length != 0)
                     break;
             }

             width       =  p->work_section[7].value;
             Aera        =  p->work_section[10].value;
             node_mass   =  density*Aera*elmt_length/8.0;
             node_Jx     =  density*elmt_length*thickness*width*width*width/12.0/8.0;
             node_Jy     =  density*Aera*elmt_length*elmt_length*elmt_length/12.0/8.0;
#ifdef DEBUG
            printf(" width      = %lf\n", width);
            printf(" length     = %lf\n", elmt_length);
            printf(" Aera       = %lf\n", Aera);
            printf(" node_mass  = %lf\n", node_mass);
            printf(" node_Jx    = %lf\n", node_Jx);
            printf(" node_Jy    = %lf\n", node_Jy);
#endif
             
             for(i = 1; i <= p->dof_per_node; i++) {
                 for (j = 1; j <= p->nodes_per_elmt; j++) {
                      k = p->dof_per_node*(j-1) + i; 
                      if(i <= 3)
                         mass->uMatrix.daa[k-1][k-1] = node_mass;
                      else{
                         if(i == 4) 
                         mass->uMatrix.daa[k-1][k-1] = node_Jx;
                         if(i == 5) 
                         mass->uMatrix.daa[k-1][k-1] = node_Jy;
                      }
                 }
             }
             
	     break;
	case CONSISTENT:

        /* MASS : [M] = [M1]+[M3] */

             e1_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
             e2_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);
             e3_ptr = MatrixAllocIndirectDouble(3,p->nodes_per_elmt);

             Mass = MatrixAllocIndirectDouble(p->size_of_stiff, p->size_of_stiff); 
             N1             = MatrixAllocIndirectDouble(3, p->size_of_stiff);
             N2             = MatrixAllocIndirectDouble(3, p->size_of_stiff);
             N1_Trans       = MatrixAllocIndirectDouble(p->size_of_stiff, 3);
             N2_Trans       = MatrixAllocIndirectDouble(p->size_of_stiff, 3);
             shp            = MatrixAllocIndirectDouble(3,8);
             x_integ_coord  = dVectorAlloc(16);
             y_integ_coord  = dVectorAlloc(16);
             x_integ_pts    = p->integ_ptr->surface_pts;
             y_integ_pts    = p->integ_ptr->surface_pts;

             z_integ_pts    = p->integ_ptr->thickness_pts;
             z_integ_coord  = dVectorAlloc(z_integ_pts + 1);
             z_weight       = dVectorAlloc(z_integ_pts + 1);


             no_integ_pts   = 0;

             gauss(z_integ_coord, z_weight,z_integ_pts);

             pgauss(x_integ_pts, &no_integ_pts, x_integ_coord, y_integ_coord, weight);

             Lamina_Sys_8node(p, e1_ptr, e2_p