elmt_psps.c

有限元分析源代码

            /* calculate body_force[] at gaussian points : body_force = sum 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*jacobain;

                    /* 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] *= jacobain;

            /* 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];
                }
            }

           /* [e] CALCULATE EQUIVALENT NODAL FORCE DUE TO SURFACE TRACTION   */

           /* add code later */

          /* 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];
                 }
          }
            
   /* ------------ EQUIVALENT NODAL LOAD UNITS ------------*/
   /* Young's Modulus E is in SI then Use SI, else use US  */
   /* ---------------------------------------------------- */
    switch(UNITS_SWITCH) {
      case ON:
       /* Initiation of Equivalent nodal load Units Buffer */

       if(UNITS_TYPE == 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 OFF:
      break;
      default:
      break;
    }
     
#ifdef DEBUG
       printf("*** in elmt_psps() : leaving case EQUIV_NODAL_LOAD: \n");
#endif
           break;

      case STRESS_UPDATE:
           break;
      case STRESS:
      case LOAD_MATRIX:
	
#ifdef DEBUG
       printf("*** in elmt_psps() : start case STRESS: or case LOAD_MATRIX: \n");
#endif

           lint = (int ) 0;
           if(isw == STRESS)      l=  no_stress_pts; /* stress pts         */
           if(isw == LOAD_MATRIX) l=  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(l*l != lint)
              pgauss(l,&lint,sg,tg,wg);

           /* element stress, strains and forces */

            if(p->no_dimen == 2 && isw == STRESS) {

                printf("\n STRESS in  Element No  %d \n", p->elmt_no);
                printf(" ================================================================================================== \n");
                printf(" Gaussion    xi        eta         x         y          Stress-11       Stress-22        Stress-12 \n");
                if(UNITS_SWITCH == OFF)
                   printf("  Points \n");
                if(UNITS_SWITCH == ON){
                   if(UNITS_TYPE == SI)
                      dp_stress = DefaultUnits("Pa");
                   else
                      dp_stress = DefaultUnits("psi");
                   printf("  Points                           %s         %s             %s             %s               %s\n", 
                             p->coord[0][0].dimen->units_name,
                             p->coord[1][0].dimen->units_name,
                             dp_stress->units_name,
                             dp_stress->units_name,
                             dp_stress->units_name);
                   free((char *) dp_stress->units_name);
                   free((char *) dp_stress);
                }
                printf(" --------------------------------------------------------------------------------------------------\n \n");
            }
#ifdef DEBUG
            if(p->no_dimen == 2 && isw == LOAD_MATRIX) {
                printf("\n NODAL LOAD in  Element No  %d \n", p->elmt_no);
                printf(" ===================================================================================\n");
                printf(" Gaussion    xi        eta         x         y          nodal_no        nodal_load  \n");
                printf("  Points \n");
                printf(" ----------------------------------------------------------------------------------- \n");
            }
#endif

           for( l= 1; l<= lint; l++) {
             /* element shape function and their derivatives */
                shape(sg[l-1], tg[l-1],p->coord,shp,&jacobain,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];
            }
            
            Mater_matrix = MATER_MAT_PLANE(E.value, nu);
            
             /* 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);

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

                  /* compute eqivalent node forces due to stress */
                  /* f_equiv = int [B]^T stress dV               */ 

                   dv = jacobain*wg[l-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);
                }
                else { /* case STRESS */

                  /* output stresses  */
                  printf("   %d  %10.4f %10.4f %10.4f %10.4f", l, sg[l-1], tg[l-1], xx, yy);
                  printf("\t%10.4f\t%10.4f\t%10.4f\n", stress[0][0], stress[1][0], stress[2][0]);  
                }
               
                for(i = 1; i <= dof_per_elmt; i++) 
                   p->nodal_loads[i-1].value += nodal_load[i-1][0];
            }

   /* ------------NODAL LOAD UNITS ------------------------*/
   /* The units type is determined by the SetUnitsType()   */
   /* ---------------------------------------------------- */
    switch(UNITS_SWITCH) {
      case ON:
        if(UNITS_TYPE == SI)
           d1 = DefaultUnits("N");
        else
           d1 = DefaultUnits("lbf");

        /* node no 1 */
        UnitsCopy( p->nodal_loads[0].dimen, d1 );
        UnitsCopy( p->nodal_loads[1].dimen, d1 );
        /* node no > 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->nodal_loads[k-1].dimen, d1 );
            }
        }
        free((char *) d1->units_name);
        free((char *) d1);
      break;
      case OFF:
      break;
      default:
      break;
    }

 /* ------------====================-------------------- */
            break;
       case MASS_MATRIX:
             
            /* compute consistent mass matrix      */
            /* p->type should be -1 for consistent */

            l = no_integ_pts; 

           printf(" no_integ_pts = %d \n", l);
           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(l*l != lint)
               pgauss(l,&lint,sg,tg,wg);

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

             /* compute shape functions */

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

                dv = density.value * wg[l-1] * jacobain;

             /* for each node compute db = shape * dv  */
                j1 = 1;
                for(j = 1; j<= p->nodes_per_elmt; j++){
                    w11 = shp[2][j-1] * dv;

                 /* compute lumped mass */
                 /* ?? store lumped mass in p->nodal_loads ?? */
                    p->nodal_loads[j1-1].value = p->nodal_loads[j1-1].value + w11; 

                 /* for each node k compute mass matrix ( upper triangular part ) */
                    k1 = j1;
                    for(k = j; k <= p->nodes_per_elmt; k++) {
                        stiff[j1-1][k1-1] += shp[2][k-1] * w11;
                        k1 = k1 + p->dof_per_node;
                    }
                    j1 = j1 + p->dof_per_node;
                } 
                      
                for (i = 1; i <= dof_per_elmt; i++) {
                   for (j = 1; j <= dof_per_elmt; j++) {
                      p->stiff->uMatrix.daa[i-1][j-1] += stiff[i-1][j-1];
                   }
                }
                 
            }

            /* compute missing parts and lower part by symmetries */

            dof_per_elmt = p->nodes_per_elmt* p->dof_per_node;

            for(j = 1; j <= dof_per_elmt; j++){
                /* ??????? Store lumped mass in p->nodal_loads ?? */
                p->nodal_loads[j].value = p->nodal_loads[j-1].value;
                for(k = j; k <= p->dof_per_node; k = k + p->dof_per_node) {
                    p->stiff->uMatrix.daa[j][k]      = p->stiff->uMatrix.daa[j-1][k-1];
                    p->stiff->uMatrix.daa[k-1][j-1]  = p->stiff->uMatrix.daa[j-1][k-1];
                    p->stiff->uMatrix.daa[k][j]      = p->stiff->uMatrix.daa[j-1][k-1];
                }  
            }

#ifdef DEBUG
      /* check element mass : sum of row of Mass */

       for (i = 1; i <= dof_per_elmt; i++)
          sum_row[i-1] =  0.0;

       for (i = 1; i <= dof_per_elmt; i++)
          for (j = 1; j <= dof_per_elmt; j++)
            sum_row[i-1] +=  p->stiff->uMatrix.daa[i-1][j-1];

       sum = 0;
       for (i = 1; i <= dof_per_elmt; i++){
         printf(" Mass M_e : sum of row[%d] = %lf \n", i, sum_row[i-1]);
         sum += sum_row[i-1];
       }
         printf(" Mass M_e : sum = %lf \n",sum);


#endif

 /* ------------ MASS UNITS ---------------------------- */
 /* Initiation of Mass Units Buffer                      */

    switch(UNITS_SWITCH) {
      case ON:
        if(UNITS_TYPE == SI || UNITS_TYPE == SI_US ) {
           d1 = DefaultUnits("Pa");
           d1 = DefaultUnits("m");
        }
        else {
           d1 = DefaultUnits("psi");
           d1 = DefaultUnits("in");
        }
        d3 = DefaultUnits("sec");

        /* node no 1 */
        UnitsMultRep( &(p->stiff->spColUnits[0]), d1, d2 );
        UnitsCopy( &(p->stiff->spColUnits[1]), &(p->stiff->spColUnits[0]) );

        UnitsPowerRep( &(p->stiff->spRowUnits[0]), d3, 2.0, NO );
        UnitsCopy( &(p->stiff->spRowUnits[1]), &(p->stiff->spRowUnits[0]) );

        /* node no > 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->stiff->spColUnits[k-1]), &(p->stiff->spColUnits[0]) );
                UnitsCopy( &(p->stiff->spRowUnits[k-1]), &(p->stiff->spRowUnits[0]) );
            }
        }
        free((char *) d1->units_name);
        free((char *) d1);
        free((char *) d2->units_name);
        free((char *) d2);
        free((char *) d3->units_name);
        free((char *) d3);
      break;
      case OFF:
      break;
      default:
      break;
    }
     
 /* ------------====================-------------------- */

#ifdef DEBUG
       printf("*** leaving elmt_psps() case : MASS_MATRIX  \n");
#endif
            break;
            default:
            break;
    }
   
    free(alpha_thermal);
    free(alpha_thermal);
    free(sum_row);
    MatrixFreeIndirectDouble(strain, 3);
    MatrixFreeIndirectDouble(stress, 3);
    MatrixFreeIndirectDouble(body_force, 3);
    MatrixFreeIndirectDouble(load, dof_per_elmt);
    MatrixFreeIndirectDouble(nodal_load, dof_per_elmt);
    MatrixFreeIndirectDouble(displ, dof_per_elmt);
    MatrixFreeIndirectDouble(stiff, dof_per_elmt);
    MatrixFreeIndirectDouble(B_matrix, 3);
    MatrixFreeIndirectDouble(B_Transpose, dof_per_elmt);
    MatrixFreeIndirectDouble(temp_matrix, dof_per_elmt);
    MatrixFreeIndirectDouble(shp, 3);
    

    return(p);
}


/* ================================================== */
/* shape function                                     */
/* ================================================== */

int shape(ss,tt,coord,shp,jacobian,no_dimen,nodes_per_elmt,node_connect,flg)
double  ss, tt, **shp,*jacobian, *node_connect;
QUANTITY                          **coord;
int              no_dimen, nodes_per_elmt, flg;
{
int i,j,k;
double  s[4], t[4], xs[2][2], sx[2][2], tp;
double  **shp_temp;
 
    shp_temp = MatrixAllocIndirectDouble(2, 4);

    s[0] =  0.5; s[1] = -0.5;
    s[2] = -0.5; s[3] =  0.5;