elmt_psps.c

利用语言编写的有限元分析软件

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

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

  switch(nodes_per_elmt) { 
    case 3: 
    case 4:

    /* form 4-node quadrilateral shape function                  */
    /* shape function:                  Ni = shape[2][i]         */
    /*                                  node no. i = 1, .. 4     */
    /* derivatives of shape functions:  dNi/d(ss) = shape[0][i]  */
    /*                                  dNi/d(tt) = shape[1][i]  */

    for(i = 1; i <= 4; i++){
        shp[2][i-1]      = (0.5 + s[i-1] * ss) * ( 0.5 + t[i-1] * tt); 
        shp_temp[0][i-1] = s[i-1] * (0.5 + t[i-1] * tt);                
        shp_temp[1][i-1] = t[i-1] * (0.5 + s[i-1] * ss);
    }

    /* form triangle by adding third and fourth node together */

    if(nodes_per_elmt == 3) { 
        shp[2][2] = shp[2][2] + shp[2][3];
        for(i = 0; i <= 1; i++)
            shp_temp[i][2] = shp_temp[i][2] + shp_temp[i][3];
    }

     /* construct jacobian matrix and its determinant */

     for(i = 1; i <= no_dimen; i++)      /* no_dimen = 2 */
         for(j = 1; j <= no_dimen; j++) {
             xs[i-1][j-1] = 0.0;
             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 Jacobain 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;


    /* form global derivatives */

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

    case 5: case 6: case 7: case 8:
       shp0(ss, tt, shp, node_connect, nodes_per_elmt);
    break;

    default:
    break;
  }
#ifdef DEBUG

   printf(" in shape () \n");
    for(j = 1; j <= nodes_per_elmt; j++){
            printf("   dN/dx[%d] = %lf ", j , shp[0][j-1]);
            printf("   dN/dy[%d] = %lf ", j , shp[1][j-1]);
            printf("       N[%d] = %lf \n", j , shp[2][j-1]);
   }
   printf(" leaving shape () \n");
#endif

   MatrixFreeIndirectDouble(shp_temp, 2);
    return;
}

/* ================================ */
/* shap0                            */
/* ================================ */

int shp0(s,t,shp,ix,nel)
double s,t,shp[4][10];
int  *ix;
int nel;
{
double s2, t2;
int i,j,k,l;

    s2 = (1- s * s)/2;
    t2 = (1 - t * t)/2;

    for(i = 5; i<= 9 ;i++)
    for(j =1;j<= 3 ; j++)
        shp[j][i] = 0.0;

    /* midside nodes (serendipty)  */

    if(ix[5] == 0) goto NEXT1;
        shp[1][5] = -s*(1-t);
        shp[2][5] = -s2;
        shp[3][5] =  s2*(1-t);

    NEXT1:
        if(nel < 6) goto NEXT7;
        if(ix[6] == 0) goto NEXT2;

        shp[1][6] = t2;
        shp[2][6] = - t *(1+s);