elmt_psps.c

有限元程序

            if(p->no_dimen == 2 && isw == STRESS) {
               if(UNITS_SWITCH == ON) {
                  free( dp_length->units_name );
                  free( dp_length );
                  free( dp_stress->units_name );
                  free( dp_stress );
               }
            }

            /* 
             *  =================================================================
             *  Nodal Load Units : units type is determined by the SetUnitsType()                 
             *  =================================================================
             */ 

            if( UNITS_SWITCH == ON ) {
                if( CheckUnitsType() == 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 STRESS_MATRIX:           /* save element stresses in working array */

            lint = (int ) 0;
            if(isw == STRESS_MATRIX)     
               ii = no_stress_pts;   /* stress 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");
               }
            }

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

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

                /* Save element level stresses in working array */
                /* Set column buffer units                      */

                if(ii == 1 ) {
                   UnitsCopy( &(p->stress->spColUnits[0]), dp_length ); 
                   UnitsCopy( &(p->stress->spColUnits[1]), dp_length ); 
                   UnitsCopy( &(p->stress->spColUnits[2]), dp_stress ); 
                   UnitsCopy( &(p->stress->spColUnits[3]), dp_stress ); 
                   UnitsCopy( &(p->stress->spColUnits[4]), dp_stress ); 
                }

                /* Zero out row buffer units */

                ZeroUnits( &(p->stress->spRowUnits[ii-1]) );

                /* Transfer xx and yy coordinates to to working stress matrix */

                p->stress->uMatrix.daa[ii-1][0] = xx;
                p->stress->uMatrix.daa[ii-1][1] = yy;

                /* Transfer internal stresses to working stress matrix */

                p->stress->uMatrix.daa[ii-1][2] = stress[0][0];
                p->stress->uMatrix.daa[ii-1][3] = stress[1][0];
                p->stress->uMatrix.daa[ii-1][4] = stress[2][0];
            }
            break;
       case MASS_MATRIX:
             
            /*
             *  ==================================================================
             *  Compute consistent mass matrix p->type should be -1 for consistent
             *  ==================================================================
             */

            /* Zero contents of integration points arrays */

            ii = no_integ_pts; 
            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;
            }

            /* Get Gauss integration points */

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

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

                /* Compute shape functions */

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

                dv = density.value * wg[ii-1] * jacobian;

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

            /* Units for Mass Matrix */

            if(UNITS_SWITCH == ON) {
               if( CheckUnitsType() == SI || CheckUnitsType() == 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;
       default:
           break;
    }

    /* [d] : free memory and leave */
   
    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() : shape function                                    
 * 
 *     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] 
 * 
 *  Input  :  double ss        -- 
 *            double tt        -- 
 *            QUANTITY **coord -- 
 *  Output : 
 *  ========================================================================= 
 */ 

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;
int nodes_per_elmt;
int            flg;
{
double  s[4], t[4], xs[2][2], sx[2][2], tp;
double  **shp_temp;
int i, j, k;

    /* [a] : Initialize arrays */
 
    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;

    /* [b] : form shape functions */

    switch(nodes_per_elmt) { 
        case 3: case 4:
           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;