elmt_frame2d.c

有限元程序

             /* Add FEF if elmt loaded */
 
             if( p->elmt_load_ptr != NULL) { 
                 p = sld07(p, STRESS);

                 /* Add FEF to joint forces */

                 fx1 = fx1  - p->nodal_loads[0].value;
                 fy1 = fy1  - p->nodal_loads[1].value;
                 mz1 = mz1  - p->nodal_loads[2].value;
                 fx2 = fx2  - p->nodal_loads[3].value;
                 fy2 = fy2  - p->nodal_loads[4].value;
                 mz2 = mz2  - p->nodal_loads[5].value;
             } 

             /* Assign force values */

             if( UNITS_SWITCH == ON ) {
                if( CheckUnitsType() == SI) {
                    dp_length = DefaultUnits("m");
                    dp_force  = DefaultUnits("N");
                }
                else if( CheckUnitsType() == US) {
                    dp_length = DefaultUnits("in");
                    dp_force  = DefaultUnits("lbf");
                }
            
                /* node no 1 */

                UnitsCopy( p->nodal_loads[0].dimen, dp_force );
                UnitsCopy( p->nodal_loads[1].dimen, dp_force );
                UnitsMultRep( p->nodal_loads[2].dimen, dp_force, dp_length );

                /* node no = 2 */

                UnitsCopy( p->nodal_loads[3].dimen, p->nodal_loads[0].dimen );
                UnitsCopy( p->nodal_loads[4].dimen, p->nodal_loads[1].dimen );
                UnitsCopy( p->nodal_loads[5].dimen, p->nodal_loads[2].dimen );
             }

             p->nodal_loads[0].value = fx1; 
             p->nodal_loads[1].value = fy1; 
             p->nodal_loads[2].value = mz1;
             p->nodal_loads[3].value = fx2; 
             p->nodal_loads[4].value = fy2; 
             p->nodal_loads[5].value = mz2; 

             /* Transfer nodal coordinates and forces/stresses to working array */
             /* Set column buffer units */

             UnitsCopy( &(p->stress->spColUnits[0]), dp_length ); 
             UnitsCopy( &(p->stress->spColUnits[1]), dp_length ); 
             UnitsCopy( &(p->stress->spColUnits[2]), dp_force ); 
             UnitsCopy( &(p->stress->spColUnits[3]), dp_force ); 
             UnitsMultRep( &(p->stress->spColUnits[4]), dp_force, dp_length );

             /* Zero out row buffer units */

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

             /* Transfer coordinates to working stress matrix */

             p->stress->uMatrix.daa[0][0] = p->coord[0][0].value;
             p->stress->uMatrix.daa[0][1] = p->coord[1][0].value;

             p->stress->uMatrix.daa[1][0] = p->coord[0][1].value;
             p->stress->uMatrix.daa[1][1] = p->coord[1][1].value;

             /* Transfer internal loads to working stress matrix */

             p->stress->uMatrix.daa[0][2] = p->nodal_loads[0].value;
             p->stress->uMatrix.daa[0][3] = p->nodal_loads[1].value;
             p->stress->uMatrix.daa[0][4] = p->nodal_loads[2].value;
             p->stress->uMatrix.daa[1][2] = p->nodal_loads[3].value;
             p->stress->uMatrix.daa[1][3] = p->nodal_loads[4].value;
             p->stress->uMatrix.daa[1][4] = p->nodal_loads[5].value;

             break;
        case STIFF: /* form element stiffness */

             cs = p->coord[0][1].value - p->coord[0][0].value;
             sn = p->coord[1][1].value - p->coord[1][0].value;
             xl = sqrt(cs * cs + sn * sn);
             cs = cs/xl; 
             sn = sn/xl;
             p->length.value = xl;

             if( UNITS_SWITCH==ON )
                 p->stiff->spColUnits[0].units_type = CheckUnitsType();

             if(p->nodes_per_elmt == 2) /* elastic elment use 2-node elmt */
                p->stiff = beamst(p, p->stiff, EA, EIzz, xl, cs, sn,
                                  p->size_of_stiff, p->dof_per_node);

             break;
        case MASS_MATRIX:

             cs = p->coord[0][1].value - p->coord[0][0].value;
             sn = p->coord[1][1].value - p->coord[1][0].value;
             xl = sqrt(cs * cs + sn * sn);
             cs = cs/xl; 
             sn = sn/xl;
             p->length.value = xl;

             /* Calculate mass = m_bar = mass/length                   */
             /* in units of (kg/m) or ((lbf-sec^2/in)/in)=(lb/in)      */
             /* if no units, then assume gravity g = 9.80665 m/sec^2   */

             if( weight != 0.0 )
               mass = weight/9.80665;
             else
               if( density.value > 0 )  mass = A * density.value ;
             else {
                  printf("ERROR >> In input file \n");
                  printf("ERROR >> You need a density value to calculate mass matrix\n");
                  exit(1);
             }

             if( UNITS_SWITCH == ON ) 
                p->stiff->spColUnits[0].units_type = CheckUnitsType();

             p->stiff = beamms(p, p->stiff, p->type, mass, xl, cs, sn,
                               p->size_of_stiff, p->dof_per_node);
             break;
        default:
             break;
    }

    return(p);
}


/* 
 *  ======================================================================
 *  beamst() : Compute beam stiffness matrix. Each element has two nodes,
 *             each having two translational and one rotation d.o.f.
 * 
 *  Input  : ARRAY  *p         -- working data structure ARRAY.
 *         : MATRIX *s         -- pointer to stiffness matrix.
 *         : double EA         -- rigidity in axial direction.
 *         : double EI         -- fexuraly rigidity.
 *         : double length     -- length of the element.
 *         : double cs         -- cosine of the element orientation.
 *         : double sn         -- sine of the element orientation.
 *         : int size_of_stiff -- size of the stiffness matrix.
 *         : int no_dof        -- 
 *  Output : MATRIX *          -- pointer to stiffness matrix.
 *  ======================================================================
 */ 

#ifdef __STDC__
MATRIX *beamst(ARRAY *p, MATRIX *s, double EA, double EI, double length,
               double cs, double sn, int size_of_stiff, int no_dof)
#else
MATRIX *beamst(p, s, EA, EI, length, cs, sn, size_of_stiff, no_dof)
ARRAY  *p;
MATRIX *s;
double  EA, EI, length, cs, sn;
int     size_of_stiff, no_dof ;
#endif
{
int             i, j, k;
int    length1, length2;
double                t;
DIMENSIONS     *d1, *d2;
DIMENSIONS     *d3, *d4;

    /* [a] : initialize working variables */

    i = no_dof + 1;
    j = no_dof + 2; 
    k = no_dof + 3;

    /* [b] : fill-in element of the stiffness matrix */

    t = EA/length; 
 
    s->uMatrix.daa[0][0]     =  t;
    s->uMatrix.daa[i-1][i-1] =  t;
    s->uMatrix.daa[0][i-1]   = -t;
    s->uMatrix.daa[i-1][0]   = -t;

    t = 12 * EI /(length*length*length);

    s->uMatrix.daa[1][1]     =  t;
    s->uMatrix.daa[j-1][j-1] =  t;
    s->uMatrix.daa[1][j-1]   = -t;
    s->uMatrix.daa[j-1][1]   = -t;

    t = (EI+ EI) / length;

    s->uMatrix.daa[2][2]   = s->uMatrix.daa[k-1][k-1] = t+t;
    s->uMatrix.daa[2][k-1] = s->uMatrix.daa[k-1][2]   = t;

    t = 6 * EI/length/length;
    s->uMatrix.daa[1][2]   = s->uMatrix.daa[2][1]       = t;
    s->uMatrix.daa[1][k-1] = s->uMatrix.daa[k-1][1]     = t;
    s->uMatrix.daa[2][j-1] = s->uMatrix.daa[j-1][2]     = -t;
    s->uMatrix.daa[j-1][k-1] = s->uMatrix.daa[k-1][j-1] = -t;

    /* [c] : fill-in the units buffer */

    if( CheckUnits() == ON ) {
       if( CheckUnitsType() == SI) {
           d1 = DefaultUnits("Pa");
           d2 = DefaultUnits("m");
           d3 = DefaultUnits("N");
           d4 = DefaultUnits("rad");
       }
       else {
           d1 = DefaultUnits("psi");
           d2 = DefaultUnits("in");
           d3 = DefaultUnits("lbf");
           d4 = DefaultUnits("rad");
       }

       /* Column buffer units */

       UnitsMultRep( &(s->spColUnits[0]), d1, d2 );
       UnitsCopy( &(s->spColUnits[1]), &(s->spColUnits[0]) );
       UnitsDivRep( &(s->spColUnits[2]), d3, d4, NO );

       UnitsCopy( &(s->spColUnits[3]), &(s->spColUnits[0]) ); 
       UnitsCopy( &(s->spColUnits[4]), &(s->spColUnits[1]) ); 
       UnitsCopy( &(s->spColUnits[5]), &(s->spColUnits[2]) ); 

       /* Row buffer units */

       ZeroUnits( &(s->spRowUnits[0]) );
       ZeroUnits( &(s->spRowUnits[1]) );
       UnitsCopy( &(s->spRowUnits[2]), d2 ); 

       UnitsCopy( &(s->spRowUnits[3]), &(s->spRowUnits[0]) ); 
       UnitsCopy( &(s->spRowUnits[4]), &(s->spRowUnits[1]) ); 
       UnitsCopy( &(s->spRowUnits[5]), &(s->spRowUnits[2]) ); 

       /* Free working memory */

       free((char *) d1->units_name);
       free((char *) d1);
       free((char *) d2->units_name);
       free((char *) d2);
       free((char *) d3->units_name);
       free((char *) d3);
       free((char *) d4->units_name);
       free((char *) d4);
    }

    /* [d] : rotate element to global coordinate system */

    s->uMatrix.daa = (double **) rotate(s->uMatrix.daa, cs, sn, size_of_stiff, no_dof);

    return(s);
}


/*
 *  ======================================================================
 *  beamms () : compute mass matrix for beam.
 *
 *  Input  : ARRAY  *p         -- working data structure ARRAY.
 *         : MATRIX *s         -- pointer to stiffness matrix.
 *         : int  mtype        -- mass matrix type (lumped or consistent).
 *         : double mass       -- lumped nodal mass.
 *         : double xl         -- length of the element.
 *         : double cs         -- cosine of the element orientation.
 *         : double sn         -- sine of the element orientation.
 *         : int nst           -- 
 *         : int ndf           -- 
 *  Output : MATRIX *          -- pointer to stiffness matrix.
 *  ======================================================================
 */

#ifdef __STDC__
MATRIX *beamms(ARRAY *p, MATRIX *s, int mtype, double mass,
               double xl, double cs, double sn, int nst, int ndf)
#else
MATRIX *beamms(p, s, mtype, mass, xl, cs, sn, nst, ndf)
ARRAY  *p;
MATRIX *s;
double mass, xl, cs, sn;  /* mass, length, cosine, ans sine */
int nst, ndf,   mtype;
#endif
{
int    i, j, k, n, i1, j1;
int      length1, length2;
double  t, s1, s2, s3, dv;
double sg[5], ag[5], ba[3], rba[3], bb[5], rbb[5];
DIMENSIONS  *d1, *d2, *d3;

    /* [a] : Assemble lumped and consistent mass matrices */

    t = mass * xl/2;
    switch(mtype) {
	case LUMPED:
             s->uMatrix.daa[0][0] = s->uMatrix.daa[1][1]  =  t;
             s->uMatrix.daa[2][2] =  t*xl*xl/12.0;   /* Using 16.0 in original version */

             s->uMatrix.daa[ndf][ndf] = t;
             s->uMatrix.daa[ndf+1][ndf+1] = t;
             s->uMatrix.daa[ndf+2][ndf+2] = t*xl*xl/12.0;   /* Using 16.0 in original version */

	     break;
	case CONSISTENT:

             gauss(sg, ag, 4);

             for(n =1; n<=4; n++) {

                 dv = 0.5 * xl * mass * ag[n];

                 s1 = (1 + sg[n])/2;
                 s2 = s1 * s1;
                 s3 = s1 * s2;

                 ba[1] = 1 - s1;
                 ba[2] = s1;
                 bb[1] = 1 - 3*s2 + s3 + s3;
                 bb[2] = xl* (s1- s2 - s2 + s3);
                 bb[3] = 3*s2 - s3 - s3;
                 bb[4] = xl  *  (-s2 + s3 );

                 rba[1] = ba[1] * dv;
                 rba[2] = ba[2] * dv;
                 for(i=1; i<= 4; i++)
                     rbb[i] = bb[i] * dv;

                 i = ndf;
                 s->uMatrix.daa[0][0] = s->uMatrix.daa[0][0] + ba[1] * rba[1];
                 s->uMatrix.daa[0][i] = s->uMatrix.daa[0][i] + ba[1] * rba[2];
                 s->uMatrix.daa[i][i] = s->uMatrix.daa[i][i] + ba[2] * rba[2];

                 for(i=1; i<= 4; i++){
                     i1 = i + 1;
                     if(i1 > 3 ) i1 = i1 + 1;
                     for(j=1; j<= 4; j++){
                         j1 = j  + 1;
                         if(j1 > 3) j1 = j1 +1;
                         s->uMatrix.daa[i1-1][j1-1] = s->uMatrix.daa[i1-1][j1-1] + bb[i]* rbb[j];
                     }
                 }
             }

             for(i = 0; i < 4; i++)
                 for(j = 0; j < i; j++)
                     s->uMatrix.daa[i][j] = s->uMatrix.daa[j][i];

	     break;
	default:
             FatalError("In elmt_frame2d() : beamms() : Type of Mass Matrix Undefined",(char*)NULL);
	     break;
    }

    /* [b] : Setup units buffers for mass matrix */

    if( CheckUnits() == ON ) {
       if( CheckUnitsType() == SI) {
           d1 = DefaultUnits("Pa");
           d2 = DefaultUnits("m");
       }
       else {
           d1 = DefaultUnits("psi");
           d2 = DefaultUnits("in");
       }
       d3 = DefaultUnits("sec");

       UnitsMultRep( &(s->spColUnits[0]), d1, d2 );
       UnitsCopy( &(s->spColUnits[1]), &(s->spColUnits[0]) );
       UnitsMultRep( &(s->spColUnits[2]), &(s->spColUnits[0]), d2 );

       UnitsPowerRep( &(s->spRowUnits[0]), d3, 2.0, NO );
       UnitsCopy( &(s->spRowUnits[1]), &(s->spRowUnits[0]) );
       UnitsMultRep( &(s->spRowUnits[2]), d2, &(s->spRowUnits[0]) );

       UnitsCopy( &(s->spColUnits[3]), &(s->spColUnits[0]) );
       UnitsCopy( &(s->spColUnits[4]), &(s->spColUnits[1]) );
       UnitsCopy( &(s->spColUnits[5]), &(s->spColUnits[2]) );

       UnitsCopy( &(s->spRowUnits[3]), &(s->spRowUnits[0]) );
       UnitsCopy( &(s->spRowUnits[4]), &(s->spRowUnits[1]) );
       UnitsCopy( &(s->spRowUnits[5]), &(s->spRowUnits[2]) );

       free((char *) d1->units_name);
       free((char *) d1);
       free((char *) d2->units_name);
       free((char *) d2);
       free((char *) d3->units_name);
       free((char *) d3);
    }

    /* [c] : rotate mass matrix to frame coordinate scheme */

    s->uMatrix.daa = (double **) rotate(s->uMatrix.daa,cs,sn,nst,ndf);

    return(s);
}