elmt_shell_8n.c

有限元程序

      for(j = 1; j <= p->nodes_per_elmt; j++) {
          k = p->dof_per_node*(j-1)-1;

          B_matrix[0][k+1] = a[j-1][0];
          B_matrix[0][k+2] = 0.0;
          B_matrix[0][k+3] = 0.0;

          B_matrix[1][k+1] = 0.0;
          B_matrix[1][k+2] = b[j-1][0];
          B_matrix[1][k+3] = 0.0;

          B_matrix[2][k+1] = 0.0;
          B_matrix[2][k+2] = 0.0;
          B_matrix[2][k+3] = c[j-1][0];

          B_matrix[3][k+1] = b[j-1][0]; 
          B_matrix[3][k+2] = a[j-1][0]; 
          B_matrix[3][k+3] = 0.0;

          B_matrix[4][k+1] = 0.0;
          B_matrix[4][k+2] = c[j-1][0];
          B_matrix[4][k+3] = b[j-1][0]; 

          B_matrix[5][k+1] = c[j-1][0];
          B_matrix[5][k+2] = 0.0;
          B_matrix[5][k+3] = a[j-1][0]; 


      /* ----------------------------------*/
      /*  Calculate Bi" matrix             */
      /*  Bi" = [ ] 6x2                    */
      /* ----------------------------------*/

          B_matrix[0][k+4] =  -d[j-1][0]*e2_ptr[0][j-1];
          B_matrix[0][k+5] =   d[j-1][0]*e1_ptr[0][j-1];
          
          B_matrix[1][k+4] =  -e[j-1][0]*e2_ptr[1][j-1];
          B_matrix[1][k+5] =   e[j-1][0]*e1_ptr[1][j-1];

          B_matrix[2][k+4] =  -g[j-1][0]*e2_ptr[2][j-1];
          B_matrix[2][k+5] =   g[j-1][0]*e1_ptr[2][j-1];

          B_matrix[3][k+4] =  -e[j-1][0]*e2_ptr[0][j-1] - d[j-1][0]*e2_ptr[1][j-1];
          B_matrix[3][k+5] =   e[j-1][0]*e1_ptr[0][j-1] + d[j-1][0]*e1_ptr[1][j-1];

          B_matrix[4][k+4] =  -g[j-1][0]*e2_ptr[1][j-1] - e[j-1][0]*e2_ptr[2][j-1];
          B_matrix[4][k+5] =   g[j-1][0]*e1_ptr[1][j-1] + e[j-1][0]*e1_ptr[2][j-1];

          B_matrix[5][k+4] =  -d[j-1][0]*e2_ptr[2][j-1] - g[j-1][0]*e2_ptr[0][j-1];
          B_matrix[5][k+5] =   d[j-1][0]*e1_ptr[2][j-1] + g[j-1][0]*e1_ptr[0][j-1];
      }

      MatrixFreeIndirectDouble(a, p->nodes_per_elmt);
      MatrixFreeIndirectDouble(b, p->nodes_per_elmt);
      MatrixFreeIndirectDouble(c, p->nodes_per_elmt);
      MatrixFreeIndirectDouble(d, p->nodes_per_elmt);
      MatrixFreeIndirectDouble(e, p->nodes_per_elmt);
      MatrixFreeIndirectDouble(g, p->nodes_per_elmt);

      MatrixFreeIndirectDouble(e1_ptr, 3);
      MatrixFreeIndirectDouble(e2_ptr, 3);
      MatrixFreeIndirectDouble(e3_ptr, 3);

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

      return(B_matrix);
}


double **Shell_Nodal_Load_8Node(nodal_load, p, co_coord, z_coord, z_integ_pt, EE, nu, task)
double                 **nodal_load;
ARRAY                            *p;
double                   **co_coord;
double                      z_coord;
int                      z_integ_pt;      /* integration point in z-direction */
double                       EE, nu;
int                            task;
{
int           i, j, k,n, ii, jj, kk;
int                       dof, size;
int                    UNITS_SWITCH;
int       surface_pts, no_integ_pts;
double               *x_integ_coord;
double               *y_integ_coord;
double         weight[16], jacobian;
double                   sum, **shp;
double            **B_matrix = NULL;
double           **B1_matrix = NULL;
double            **B1_Trans = NULL;
double           **J_inverse = NULL;
double                  **TT = NULL;
double                     **stress; 
double            **nodal_load_temp;


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

    UNITS_SWITCH = CheckUnits();

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

    TT               = MatrixAllocIndirectDouble(6, 6);
    J_inverse        = MatrixAllocIndirectDouble(3, 3);
    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);

    nodal_load_temp  = MatrixAllocIndirectDouble(p->size_of_stiff, 1);
    stress           = MatrixAllocIndirectDouble(p->dof_per_node, 1);

    size           = p->size_of_stiff;
    dof            = p->dof_per_node;
    surface_pts    = p->integ_ptr->surface_pts;
    no_integ_pts   = 0;
    pgauss(surface_pts, &no_integ_pts, x_integ_coord, y_integ_coord, weight);
   
    for(i = 1; i <= size; i++) {
      nodal_load[i-1][0]      = 0.0;
      nodal_load_temp[i-1][0] = 0.0;
    }

       if(task == STRESS && z_integ_pt == 1){ /* print elmement stress */
          if(p->elmt_no == 1)
             printf(" Element : %s \n Material : %s \n\n", p->elmt_type, p->material_name);

          printf("\n STRESS in  Element No  %d \n",p->elmt_no);
          printf(" =============================================================================================================== \n");
          printf(" Gaussion    xi       eta        gamma    Stre-xx         Stre-yy         Stre-xy         Stre-yz        Stre-zx \n");
          if(UNITS_SWITCH == OFF)
             printf("  Points \n");
       }


    for( ii = 1; ii <= no_integ_pts; ii++) {

       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 <= p->size_of_stiff; j++) {
               B1_matrix[i-1][j-1] = 0.0;
               B1_Trans[j-1][i-1]  = 0.0;
               for(k = 1; k <= p->dof_per_node + 1; k++)
                   B1_matrix[i-1][j-1] += TT[i-1][k-1]*B_matrix[k-1][j-1];
           }

       for(i = 3; i <= p->dof_per_node; i++)
           for(j = 1; j <= p->size_of_stiff; j++) {
               B1_matrix[i-1][j-1] = 0.0;
               B1_Trans[j-1][i-1]  = 0.0;
               for(k = 1; k <= p->dof_per_node + 1; k++)
                   B1_matrix[i-1][j-1] += TT[i][k-1]*B_matrix[k-1][j-1];

           }

       for( i = 1; i <= p->dof_per_node; i++) 
           for(j = 1; j <= p->size_of_stiff; j++) 
               B1_Trans[j-1][i-1] = B1_matrix[i-1][j-1];

       /* update the stress in array pointer */

       kk = no_integ_pts*(z_integ_pt-1)+ii;

#ifdef DEBUG
       printf("jj = %d, kk = %d z_integ_pt = %d, no_integ_pts = %d \n",
               ii, kk, z_integ_pt, no_integ_pts);
       printf(" In Shell_Nodal_Load_8Node(): begin Stress_Update_8Node(): \n");
#endif
    
       Stress_Update_8Node(p, co_coord, B1_matrix, kk);

#ifdef DEBUG
       printf(" end of Stress_Update_8Node() \n");
#endif

       /* IF THE TASK IS TO DO STRESS UPDATE or PRINT STRESS, */
       /* STOP HERE AND CONTINUE FOR THE NEXT LOOP            */
      
       if(task == STRESS_UPDATE)
          goto STRESS_UPDATE_END; 

       if(task == STRESS){ /* print elmement stress : local */

          if(UNITS_SWITCH == ON) {
             if(kk == 1) {
                printf("  Points                                    %s             %s             %s             %s            %s \n",
                     p->stress->spRowUnits[0].units_name,
                     p->stress->spRowUnits[1].units_name,
                     p->stress->spRowUnits[2].units_name,
                     p->stress->spRowUnits[3].units_name,
                     p->stress->spRowUnits[4].units_name);
                printf(" ---------------------------------------------------------------------------------------------------------------\n \n");
             }
             printf("   %d  %10.4f %10.4f %10.4f", kk, x_integ_coord[ii-1], y_integ_coord[ii-1], z_coord);
             printf("\t%11.4e\t%11.4e\t%11.4e\t%11.4e\t%11.4e\n",
                    p->stress->uMatrix.daa[0][kk-1]/p->stress->spRowUnits[0].scale_factor,
                    p->stress->uMatrix.daa[1][kk-1]/p->stress->spRowUnits[1].scale_factor,
                    p->stress->uMatrix.daa[2][kk-1]/p->stress->spRowUnits[2].scale_factor,
                    p->stress->uMatrix.daa[3][kk-1]/p->stress->spRowUnits[3].scale_factor,
                    p->stress->uMatrix.daa[4][kk-1]/p->stress->spRowUnits[4].scale_factor);
          }
          if(UNITS_SWITCH == OFF) {
             if(kk == 1) 
                printf(" ---------------------------------------------------------------------------------------------------------------\n \n");
             printf("   %d  %10.4f %10.4f %10.4f", kk, x_integ_coord[ii-1], y_integ_coord[ii-1], z_coord);
             printf("\t%11.4e\t%11.4e\t%11.4e\t%11.4e\t%11.4e\n",
                    p->stress->uMatrix.daa[0][kk-1],
                    p->stress->uMatrix.daa[1][kk-1],
                    p->stress->uMatrix.daa[2][kk-1],
                    p->stress->uMatrix.daa[3][kk-1],
                    p->stress->uMatrix.daa[4][kk-1]);
          }
          goto STRESS_UPDATE_END; 
       }

       /* calculate the element nodal_load   */

        for( i = 1; i <= p->dof_per_node; i++) {
            stress[i-1][0] = p->stress->uMatrix.daa[i-1][kk-1];
        }

        nodal_load_temp = dMatrixMultRep(nodal_load_temp, B1_Trans, size, dof, stress, dof, 1);
        
        for(i = 1; i <= p->size_of_stiff; i++) {
            nodal_load[i-1][0] += nodal_load_temp[i-1][0]*jacobian*weight[ii-1]; 
        } 

STRESS_UPDATE_END: 
        ;

    }  /* end of gaussian integration */

    free((char *) x_integ_coord);
    free((char *) y_integ_coord);
    MatrixFreeIndirectDouble(shp, 3);

    MatrixFreeIndirectDouble(nodal_load_temp, p->size_of_stiff);
    MatrixFreeIndirectDouble(B_matrix, p->dof_per_node + 1);
    MatrixFreeIndirectDouble(B1_matrix, p->dof_per_node);
    MatrixFreeIndirectDouble(B1_Trans, p->size_of_stiff);
    MatrixFreeIndirectDouble(stress, p->dof_per_node);

#ifdef DEBUG 
     dMatrixPrint("nodal_load surface ", nodal_load, p->size_of_stiff, 1);
     printf(" Leaving Shell_Nodal_Load_8Node() \n");
#endif

    return(nodal_load);
}


#ifdef __STDC__
void Lamina_Sys_8node(ARRAY *p, double **e1_ptr, double **e2_ptr, double **e3_ptr)
#else
void Lamina_Sys_8node(p, e1_ptr, e2_ptr, e3_ptr)
ARRAY                         *p;
double                  **e1_ptr;
double                  **e2_ptr;
double                  **e3_ptr;
#endif
{
double dof, size, temp;

static double x75, x86;
static double y75, y86;
static double z75, z86;

double s_norm;

static double **r75_ptr = NULL, **r86_ptr = NULL;
static double **ez_ptr  = NULL, **ey_ptr  = NULL;

static double **v1_ptr  = NULL;
static double **v2_ptr  = NULL;
static double **v3_ptr  = NULL;

double **s_ptr = NULL;

int i, j, k, ii, jj, kk, n, n1, n2, k1, k2;
    
#ifdef DEBUG
    printf(" Enter Lamina_Sys_8node() \n");
#endif

    /* =================================================== */
    /* In this case, the normal direction in each node     */
    /* is simplified as the normal direction to element    */
    /* i.e. each element is considered as a plane element  */
    /* =================================================== */

    dof  = p->dof_per_node;
    size = p->size_of_stiff;

    ez_ptr  = MatrixAllocIndirectDouble(3, 1);
    ey_ptr  = MatrixAllocIndirectDouble(3, 1);

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

    ez_ptr[0][0] = 0.0;
    ez_ptr[1][0] = 0.0;
    ez_ptr[2][0] = 1.0;

    ey_ptr[0][0] = 0.0;
    ey_ptr[1][0] = 1.0;
    ey_ptr[2][0] = 0.0;

 /* --------------------------------------------------- */
 /* [1] Orientatuions of local base vectors             */ 
 /* --------------------------------------------------- */

 /* [a] Compute e3_ptr                      */
 /*     pointer normal to the shell surface */
     
#ifdef DEBUG
    printf(" In elmt_shell_Belytschko_explicit(): \n");
    printf("    : computing e3_ptr \n");
#endif

    x75 =  p->coord[0][6].value - p->coord[0][4].value;
    y75 =  p->coord[1][6].value - p->coord[1][4].value;
    z75 =  p->coord[2][6].value - p->coord[2][4].value;

    r75_ptr = MatrixAllocIndirectDouble(3, 1);
    r75_ptr[0][0] = x75;
    r75_ptr[1][0] = y75;
    r75_ptr[2][0] = z75;

    x86 =  p->coord[0][7].value - p->coord[0][5].value;
    y86 =  p->coord[1][7].value - p->coord[1][5].value;
    z86 =  p->coord[2][7].value - p->coord[2][5].value;

    r86_ptr = MatrixAllocIndirectDouble(3, 1);
    r86_ptr[0][0] = x86;
    r86_ptr[1][0] = y86;
    r86_ptr[2][0] = z86;

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

    s_ptr   = MatrixAllocIndirectDouble(3, 1);
    s_ptr   = dVmatrixCrossProduct(s_ptr, r75_ptr, 3, 1, r86_ptr, 3,1);

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

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