elmt_frame3d.c

有限元程序

	            if(n > m)                        /* sym terms */
	               s->uMatrix.daa[m-1][n-1] = s->uMatrix.daa[n-1][m-1] ;         
	            s->uMatrix.daa[n+no_dof-1][m+no_dof-1] = s->uMatrix.daa[n-1][m-1]; /* Mkk */         
	        }
	    }

	    /* ------------------------------------ */
	    /* off diagonal terms & sign adjustment */
	    /* ------------------------------------ */

	    s->uMatrix.daa[4+no_dof][2+no_dof] = s->uMatrix.daa[2+no_dof][4+no_dof] =   t* 22 * xl;
	    s->uMatrix.daa[5+no_dof][1+no_dof] = s->uMatrix.daa[1+no_dof][5+no_dof] = - t* 22 * xl;

	    /*-------------------------------------------*/
	    /* Mkj, Mjk terms - off main diagonal terms  */
	    /*-------------------------------------------*/

            s->uMatrix.daa[no_dof][0]                          =  t* 70;
            s->uMatrix.daa[no_dof+ 1][1] = s->uMatrix.daa[no_dof+2][2]   =  t* 54;
            s->uMatrix.daa[no_dof+ 3][3]                       =  t*  70 * rg * rg;
            s->uMatrix.daa[no_dof+ 4][4] = s->uMatrix.daa[no_dof+ 5][5]  = - t* 3 * xl * xl;
            s->uMatrix.daa[no_dof+ 4][2] = s->uMatrix.daa[no_dof+ 1][5]  =   t* 13 * xl;
            s->uMatrix.daa[no_dof+ 5][1] = s->uMatrix.daa[no_dof+ 2][4]  = - t* 13 * xl;

            /* Symmmetric terms */

            for(m = 1; m <= no_dof; m++)
                for(n = 1; n <= no_dof; n++)
                    s->uMatrix.daa[n-1][m+no_dof-1] = s->uMatrix.daa[n+no_dof-1][m-1];

	    break;

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


#ifdef DEBUG
    printf("elmt mass[12][12]");
    MatrixPrintIndirectDouble(s);
#endif

    /* ---------------------- */
    /* initial unit of 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]) );
       UnitsCopy( &(s->spColUnits[2]), &(s->spColUnits[0]) );
       UnitsMultRep( &(s->spColUnits[3]), &(s->spColUnits[0]), d2 );
       UnitsCopy( &(s->spColUnits[4]), &(s->spColUnits[3]) );
       UnitsCopy( &(s->spColUnits[5]), &(s->spColUnits[3]) );

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

       UnitsCopy( &(s->spColUnits[6]), &(s->spColUnits[0]) );
       UnitsCopy( &(s->spColUnits[7]), &(s->spColUnits[1]) );
       UnitsCopy( &(s->spColUnits[8]), &(s->spColUnits[2]) );
       UnitsCopy( &(s->spColUnits[9]), &(s->spColUnits[3]) );
       UnitsCopy( &(s->spColUnits[10]), &(s->spColUnits[4]) );
       UnitsCopy( &(s->spColUnits[11]), &(s->spColUnits[5]) );

       UnitsCopy( &(s->spRowUnits[6]), &(s->spRowUnits[0]) );
       UnitsCopy( &(s->spRowUnits[7]), &(s->spRowUnits[1]) );
       UnitsCopy( &(s->spRowUnits[8]), &(s->spRowUnits[2]) );
       UnitsCopy( &(s->spRowUnits[9]), &(s->spRowUnits[3]) );
       UnitsCopy( &(s->spRowUnits[10]), &(s->spRowUnits[4]) );
       UnitsCopy( &(s->spRowUnits[11]), &(s->spRowUnits[5]) );

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

    /* ---------------------- */
    /* rotate to global axis  */
    /* ---------------------- */

    s->uMatrix.daa = (double **) rotate3d(s->uMatrix.daa, rot, no_dof);

    return(s);
}


/* ====================================================== */
/* function   tmat(rot, type, p)                          */
/*    Transformation matrix for 2,3d element              */
/*    Generates a transformation matrix                   */ 
/* Input:                                                 */
/*        rot : rotation matrix (12 x 12)                 */
/*        type: elmt type                                 */
/*           p: ARRAY * storing elmt info                 */
/* ====================================================== */

#ifdef __STDC__
double **tmat(double **rot, int type, ARRAY *p)
#else
double **tmat(rot, type, p)
double **rot;
int    type;
ARRAY  *p;
#endif
{
int    i,j,n,sr,ii,jj,kk,nf, ndff;
double t[4][4];
double alpha, cx,cy,cz,csa,sna,elone, el;

    alpha = p->eangle.value;
    sr    = p->size_of_stiff;

    for(i = 0; i <= 2; i++)
        for(j = 0; j <= 2; j++)
            t[i][j] = 0.0;

    for(i = 1; i <= sr; i++)
        for(n =1; n<=sr;n++)
            rot[i-1][n-1] = 0.0;

    cx = p->coord[0][1].value - p->coord[0][0].value;
    cy = p->coord[1][1].value - p->coord[1][0].value;

    switch(type) {
        case 1: /* plane truss */
	case 2:
	case 3:
	case 4: /* plane frame, with axial force */
             el = sqrt(cx * cx + cy * cy );
             cx = cx/el; 
             cy = cy/el;
             if(type == 4) {
                t[0][0] = cx;
                t[0][1] = cy;
                t[1][0] = -cy;
                t[1][1] = cx;
                if(type == 1) break;
                   t[2][2] = 1;
             }
             else {
                t[0][0] = cx;
                t[0][2] = cy;
                t[1][1] = 1;
                t[2][0] = -cy;
                t[2][2] = cx;
             } 
             break;
        case 5: /* space truss */
        case 6: /* space frame */
             cz = p->coord[2][1].value - p->coord[2][0].value;
             el = sqrt(cx * cx + cy * cy + cz * cz);
             cx = cx/el; 
             cy = cy/el;
             cz = cz/el;
             if(type == 5){
                t[0][0] = cx*cx;
                t[0][1] = cx*cy;
                t[0][2] = cx*cz;
                t[1][1] = cy*cy;
                t[1][2] = cy*cz;
                t[2][2] = cz*cz;
                for(i = 2; i <= 3; i++){
                    ii = i -1;
                    for(j = 1; j <= ii; j++) 
                        t[i-1][j-1] = t[j-1][i-1];
                }
             }
             else {
                csa = cos(alpha);
                sna = sin(alpha);
                t[0][1] = cy;
                if(cx == 0.0 && cz == 0.0){  /* vertical space frame */
                   t[1][0] = -cy*csa;
                   t[1][2] = sna;
                   t[2][0] = cy * sna;
                   t[2][2] = csa;
                }
		else {
                   elone = sqrt(cx * cx + cz * cz);
                   t[0][0] = cx;
                   t[0][2] = cz;
                   t[1][0] = (-cx*cy*csa - cz*sna)/elone;
                   t[1][2] = (-cy*cz*csa + cx*sna)/elone;
                   t[1][1] = csa * elone;
                   t[2][0] = ( cx*cy*sna - cz*csa)/elone;
                   t[2][1] = -sna * elone;
                   t[2][2] = ( cy*cz*sna + cx*csa)/elone;
               }
            }
            break;
        default:
            break;
    }

    if(type == 6){
       nf = 4;
       ndff = 3;
    }
    else{
       nf   = 2;
       ndff = p->dof_per_node;
    }

    for(i=1;i<=nf;i++) {
        jj = (i-1) * ndff; 
        kk = (i-1) * ndff;
        for(j=1;j<=ndff;j++){
            jj = jj+ 1;
            for(n=1;n<=ndff;n++){
                kk = kk  + 1;
                rot[jj-1][kk-1] = t[j-1][n-1];
            }
            kk = (i-1) * ndff;
        }
    }

    return(rot);
}


/* ===================================================================== */
/* MATRIX *rotate3d(s, r ,ndf)                                           */
/* Rotation matrix  for 3d frame element                                 */
/* Multiplies  the stiffness matrix "s" by the transformation matrix "r" */
/* (ie. Kg = Tt*Kl*T)  Return Kg                                         */
/* ===================================================================== */

#ifdef __STDC__
double **rotate3d(double **s, double **r , int no_dof)
#else
double **rotate3d(s, r ,no_dof)
double **s, **r;
int    no_dof;
#endif
{
int    i,j,n;
double t;
double **rt, **asj, **sj, **sij, **kj;

     if(r[0][0] == 1.0) return(s);

  /* transpose of rotation matrix  r */

     rt = MatrixAllocIndirectDouble(no_dof, no_dof);
     for(i = 1; i <= no_dof; i++)
       for(n =1; n<= no_dof; n++)
           rt[n-1][i-1] = r[i-1][n-1];

  /* sii terms of transformed stiffness matrix s */

     kj = MatrixAllocIndirectDouble(no_dof, no_dof);
     for(i = 1; i <= no_dof; i++)
         for(n =1; n<= no_dof; n++)
             kj[i-1][n-1] = s[i-1][n-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<=no_dof; n++)
             s[i-1][n-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);

  /* sjj terms of transformed stiffness matrix s */

     for(i = 1; i <= no_dof; i++)
         for(n =1; n<=no_dof;n++)
             kj[i-1][n-1] = s[i+no_dof-1][n+no_dof-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<=no_dof; n++)
             s[i+no_dof-1][n+no_dof-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);

  /* sij terms of transformed stiffness matrix s */

     for(i = 1; i <= no_dof; i++)
         for(n =1; n<= no_dof; n++)
             kj[i-1][n-1] = s[i-1][n+no_dof-1];

     asj = (double **) dMatrixMult( kj, no_dof, no_dof, r, no_dof, no_dof);
     sj  = (double **) dMatrixMult( rt, no_dof, no_dof, asj, no_dof, no_dof);

     for(i=1;i<= no_dof; i++)
         for(n=1;n<= no_dof ;n++)
             s[i-1][n+no_dof-1] = s[n+no_dof-1][i-1] = sj[i-1][n-1];

     MatrixFreeIndirectDouble(asj, no_dof);
     MatrixFreeIndirectDouble( sj, no_dof);
     MatrixFreeIndirectDouble( kj, no_dof);
     MatrixFreeIndirectDouble( rt, no_dof);

     return(s);
}


/* 
 *  =============================================================
 *  print_property_frame_3d() : Print FRAME_3D Element Properties
 *  =============================================================
 */ 

#ifdef __STDC__
void print_property_frame_3d(EFRAME *frp, int i)
#else
void print_property_frame_3d(frp, i)
EFRAME    *frp;
int          i;                 /* elmt_attr_no */
#endif
{
int     UNITS_SWITCH;
ELEMENT_ATTR    *eap;

     UNITS_SWITCH = CheckUnits();
     eap = &frp->eattr[i-1];