elmt_psps.c

有限元程序

            
                /* Multiply Jacobian determinant with weight  */
                /* coefficents and material matrix            */

                jacobian *= wg[ii-1];
           
                /* [a] Calculate equivalent nodal forces due to initial strains */

                if(p->nodal_init_strain != NULL) {

                   for( i = 1; i <= 3; i++ ) 
                   for( j = 1; j <= 3; j++ )  
                        Mater_matrix[i-1][j-1] *= jacobian;

                   /* Transpose [B] matrix */

                   for(i = 1; i <= 3; i++) 
                   for(j = 1; j <= dof_per_elmt; j++)
                        B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];

                   /* Calculate strain[] at gaussian points : strain = sum Ni*nodal_strain */

                   for (i = 1; i <= 3; i++)
                      strain[i-1][0]  =  0.0;
                
                   for (i = 1; i <= 3; i++) {
                   for( j = 1; j <= p->nodes_per_elmt; j++) {
                      strain[i-1][0]  += shp[2][j-1] * p->nodal_init_strain[i-1][j-1];
                   }
                   }

                   /* [Mater]_3x3 * [strain]_3x1 and save in [stress]_3x1    */

                   stress = dMatrixMultRep(stress, Mater_matrix, 3, 3, strain, 3, 1);
           
                   /* Mutiply [B]^T * [Mater]*[strain] */

                   if(nodal_load == NULL ) { 
                      nodal_load = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                   } else {
                      load = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                      for (i = 1; i<= dof_per_elmt; i++)  {
                          nodal_load[i-1][0] += load[i-1][0];
                      }
                   }
                }

                /* [b] Calculate equivalent nodal force due to internal stress */

                if(p->nodal_init_stress != NULL) {

                   /* Transpose [B] matrix */

                   for(i = 1; i <= 3; i++)
                   for(j = 1; j <= dof_per_elmt; j++)
                       B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];

                   /* Calculate stress[] at gaussian points : stress = sum Ni*nodal_stress */
    
                   for (i = 1; i <= 3; i++)
                      stress[i-1][0]  =  0.0;
                
                   for(i = 1; i <= 3; i++) {
                      for(j = 1; j <= p->nodes_per_elmt; j++) {
                          stress[i-1][0]  += shp[2][j-1] * p->nodal_init_stress[i-1][j-1].value;
                      }
                      stress[i-1][0] *= jacobian;
                   }

                   /* Multiply [B]^T * [stress] */
               
                   if(nodal_load == NULL) {
                      nodal_load  = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                   } else {
                      load        = dMatrixMultRep(load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                      for (i = 1; i<= dof_per_elmt; i++) 
                          nodal_load[i-1][0] -= load[i-1][0];
                   }
                }

                /* [c] Calculate equivalent nodal force due to body force */

                if(p->nodal_body_force != NULL) {

                   /* Calculate body_force[] at gaussian points --  */
                   /* body_force = sum of [ Ni*body_force ]         */

                   for (i = 1; i <= p->no_dimen; i++)
                      body_force[i-1][0]  =  0.0;
                
                   for(i = 1; i <= p->no_dimen; i++) {
                      for( j = 1; j <= p->nodes_per_elmt; j++) {
                          body_force[i-1][0] 
                               += shp[2][j-1] * p->nodal_body_force[i-1][j-1].value*jacobian;

                          /* Mutiply [N]^T * [body_force] */

                          k = (j-1)*p->no_dimen + i;
                          if(nodal_load == NULL) {
                             nodal_load[k-1][0] =  shp[2][j-1]*body_force[i-1][0];
                          } else
                             nodal_load[k-1][0] += shp[2][j-1]*body_force[i-1][0];
                      }
                   }
                }

                /* [d] Calculate equivalent nodal force due to temperature change */ 
  
                if(p->nodal_init_strain != NULL) {

                   for(i = 1; i <= 3 ; i++ ) 
                   for(j = 1; j <= 3; j++)  
                       Mater_matrix[i-1][j-1] *= jacobian;

                   /* Transpose [B] matrix */

                   for(i = 1; i <= 3; i++)
                   for(j = 1; j <= dof_per_elmt; j++)
                       B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];

                   /*  [B]^T * [Mater]  and save in [B]^T    */

                   temp_matrix = dMatrixMultRep(temp_matrix, B_Transpose,
                                                dof_per_elmt, 3, Mater_matrix, 3, 3);
               
                   /* Calculate strain[] at gaussian points : */
                   /* strain = sum Ni*nodal_strain            */

                   for(i = 1; i <= 3; i++)
                       strain[i-1][0]  =  0.0;
                
                   for(j = 1; j <= p->nodes_per_elmt; j++) {
                       temperature  += shp[2][j-1] * p->nodal_temp[j-1].value;
                   }

                   for(i = 1; i <= 2; i++) {
                       strain[i-1][0]  = temperature *alpha_thermal[i-1];
                   }

                   strain[2][0] = 0.0;
                    
                   /* mutiply [B]^T * [Mater]* [strain] */
               
                   if(nodal_load == NULL) {
                      nodal_load  = dMatrixMultRep(nodal_load, temp_matrix,
                                                   dof_per_elmt, 3, strain, 3, 1);
                   } else {
                      load        = dMatrixMultRep(load, temp_matrix, dof_per_elmt, 3, strain, 3, 1);
                      for(i = 1; i<= dof_per_elmt; i++) 
                          nodal_load[i-1][0] += load[i-1][0];
                   }
                }

                /* Transfer nodal_load to p->equiv_nodal_load */

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

            /*
             *  =================================================== 
             *  Initiation of Equivalent nodal load Units Buffer   
             *  
             *  Note : If Young's Modulus E is in SI then Use SI, 
             *         otherwise, use US. 
             *  =================================================== 
             */

            if (UNITS_SWITCH == ON) {
               if( CheckUnitsType() == SI)
                   d1 = DefaultUnits("N");
               else
                   d1 = DefaultUnits("lbf");

               /* node 1 */

               UnitsCopy( &(p->equiv_nodal_load->spRowUnits[0]), d1 ); 
               UnitsCopy( &(p->equiv_nodal_load->spRowUnits[1]), d1 );

               /* node i  i > 1*/
               for(i = 2; i <= p->nodes_per_elmt; i++) {
                  for(j = 1; j <= p->dof_per_node; j++) {
                      k  = p->dof_per_node*(i-1) + j;
                      UnitsCopy( &(p->equiv_nodal_load->spRowUnits[k-1]), d1 ); 
                  }
               }

               ZeroUnits( &(p->equiv_nodal_load->spColUnits[0]) );

               free((char *) d1->units_name);
               free((char *) d1);
            }
            break;
       case STRESS_UPDATE:
            break;
       case LOAD_MATRIX:
       case STRESS:          /* compute and print element stresses */

            lint = (int ) 0;
            if(isw == STRESS)     
               ii = no_stress_pts;   /* stress pts         */
            if(isw == LOAD_MATRIX)
               ii = no_integ_pts;    /* guassian integ pts */

            /* Initilize nodal_load */
      
            for(i = 1; i <= dof_per_elmt; i++) {
               nodal_load[i-1][0] = 0.0;
               load[i-1][0]       = 0.0;
            }

            for (i = 1; i<= no_integ_pts*no_integ_pts; i++) {
               sg[i-1] = 0.0; tg[i-1] = 0.0; wg[i-1] = 0.0;
            }

            if(ii*ii != lint)
               pgauss( ii, &lint, sg, tg, wg);

            /* Get units */

            if( UNITS_SWITCH == ON ) {
               if( CheckUnitsType() == SI) {
                   dp_length = DefaultUnits("m");
                   dp_stress = DefaultUnits("Pa");
               } else {
                   dp_length = DefaultUnits("in");
                   dp_stress = DefaultUnits("psi");
               }
            }

            /* Print Element Stresses, Strains and Forces */

            if( isw == STRESS ) {

               printf("\n");
               printf("Stresses in Element No %d\n", p->elmt_no);
               printf("=======================================================================\n");
               printf("Gaussion         x          y     stress-11     stress-22     stress-12 \n");

               if(UNITS_SWITCH == ON) {
                  printf("  Points       %3s        %3s", dp_length->units_name,
                                                          dp_length->units_name );
                  printf("    %10s    %10s    %10s\n", dp_stress->units_name,
                                                       dp_stress->units_name,
                                                       dp_stress->units_name );
               } else {
                  printf("  Points \n");
               }

               printf("=======================================================================\n");
            }

            /* Compute (3x3) Constituitive Material Matrix */

            Mater_matrix = MaterialMatrixPlane( E.value, nu, dFlag );

            /* Gauss Integration */

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

                /* Get element shape function and their derivatives */

                shape( sg[ii-1], tg[ii-1], p->coord,shp, &jacobian, p->no_dimen,
                       p->nodes_per_elmt, p->node_connect,0 );

                /* Form [B] matrix */
            
                for(j = 1; j <= p->nodes_per_elmt; j++) { 
                    k = 2*(j-1);
                    B_matrix[0][k]   = shp[0][j-1];
                    B_matrix[1][k+1] = shp[1][j-1];
                    B_matrix[2][k]   = shp[1][j-1];
                    B_matrix[2][k+1] = shp[0][j-1];
                }
            
                /* Calculate strains at guassian integretion pts */

                for(i = 1;i <= 3; i++)
                    strain[i-1][0] = 0.0;

                xx = 0.0; yy = 0.0;
                for(j = 1; j <= p->nodes_per_elmt; j++) {
                    xx = xx + shp[2][j-1] * p->coord[0][j-1].value;
                    yy = yy + shp[2][j-1] * p->coord[1][j-1].value;

                    /*  converting p->displ into a array */
                    
                    for ( k = 1; k <= p->dof_per_node; k++) {
                       j1 = p->dof_per_node*(j-1) + k; 
                       displ[j1-1][0] = p->displ->uMatrix.daa[k-1][j-1];
                    }
                }

                strain = dMatrixMultRep(strain, B_matrix, 3, dof_per_elmt, displ, dof_per_elmt, 1);
                
                /* Compute Stress */

                stress = dMatrixMultRep(stress, Mater_matrix, 3, 3, strain, 3, 1);

                /* Compute equivalent nodal forces for stresses */
                /* F_equiv = integral of [B]^T stress dV        */ 

                if(isw == LOAD_MATRIX) {       
                   for(i = 1; i <= 3; i++)
                       for(j = 1; j <= dof_per_elmt; j++)
                           B_Transpose[j-1][i-1] = B_matrix[i-1][j-1];

                  dv = jacobian*wg[ii-1];

                  for (i = 1; i<= 3; i++)
                       stress[i-1][0] *= dv;

                  nodal_load = dMatrixMultRep(nodal_load, B_Transpose, dof_per_elmt, 3, stress, 3, 1);
                }

                for(i = 1; i <= dof_per_elmt; i++) 
                    p->nodal_loads[i-1].value += nodal_load[i-1][0];

                /* Print stresses  */

                if(isw == STRESS && UNITS_SWITCH == ON) {
                   printf(" %7d %10.4f %10.4f", ii,
                            xx/dp_length->scale_factor,
                            yy/dp_length->scale_factor);
                   printf(" %12.4e  %12.4e  %12.4e\n",
                            stress[0][0]/dp_stress->scale_factor,
                            stress[1][0]/dp_stress->scale_factor,
                            stress[2][0]/dp_stress->scale_factor );
                }

                if(isw == STRESS && UNITS_SWITCH == OFF) {
                   printf(" %7d %10.4f %10.4f", ii, xx, yy);
                   printf(" %12.4e  %12.4e  %12.4e\n",
                            stress[0][0], stress[1][0], stress[2][0] );
                }
            }

            /* Free memory for units used in printing stresses */