fe_matrix.c

矩阵运算〉〉 〉〉 〉〉 〉〉 。

1066: /* Action:  Ap = Tpj . Aj; Disp:    Dj = Ttpj . Dp       */
1067: /*   -- for framed structures and conn rigid bodies      */
1068: /* function double **Transformation_Matrix               */
1069: /*                    (t,type,no_dof,cx,cy,cz)              */
1070: /*  input;  t         :allocated matrix                  */
1071: /*          type      : type of element                  */
1072: /*          no_dof       : no of dof                        */
1073: /*          cx,cy,cz  : coordinates of  point j wrt      */
1074: /*                                    reference point p  */
1075: /* =======================================================*/
1076: 
1077: #ifdef __STDC__
1078: double **Transformation_Matrix(double **t, char *type, int no_dof, double cx, double cy, double cz)
1079: #else
1080: double **Transformation_Matrix(t, type,no_dof,cx,cy,cz)
1081: double **t;
1082: char   *type;
1083: int    no_dof;
1084: double cx,cy,cz;
1085: #endif
1086: {
1087: int i,j,n;
1088:      for(n = 0; n <= NO_ELEMENTS_IN_LIBRARY; n++) { 
1089: 
1090:          if(!strcmp(type, elmt_library[n].name))
1091:          break;
1092:        }
1093: 
1094:        for(i = 1; i <= elmt_library[n].no_dof; i++)
1095:        for(j = 1; j <= elmt_library[n].no_dof; j++)
1096:            t[i][j] = 0.0;
1097: 
1098:        if(elmt_library[n].no_dof == 3) {
1099:         /*2d beam:2d bar:2d frame: grid*/
1100: 
1101:              t[1][1] = 1.0;
1102:              t[2][2] = 1.0;
1103:              if(!strcmp(type,"BEAM_2D"))
1104:                 t[2][1] = cx;
1105:              else if(!strcmp(type,"TRUSS_2D")) {
1106:                 t[3][1] = -cy;
1107:                 t[3][2] =  cx;
1108:              }
1109:              else if(!strcmp(type,"FRAME_2D")){
1110:                 t[3][1] = -cy;
1111:                 t[3][2] =  cx;
1112:                 t[3][3] =  1;
1113:              }
1114:              else {                               /*  GRID_XY */
1115:                 t[1][3] =  cy;
1116:                 t[2][3] = -cx;
1117:                 t[3][3] =  1;
1118:              }
1119: 
1120:           }
1121: 
1122:        if(elmt_library[n].no_dof == 6) {
1123: 
1124:              for(i = 1; i <= 6; i++)
1125:                  t[i][i] = 1.0;
1126:      
1127:              if(!strcmp(type,"TABLE_XZ")) {   /* 3d space table o== r storey */
1128:                 t[3][1] =  cz;
1129:                 t[3][2] = -cx;
1130:              }
1131:              else if(!strcmp(type,"TRUSS_3D")) {   /* space truss */
1132:                 t[5][1] =  cz;
1133:                 t[6][1] = -cy;
1134:                 t[4][2] = -cz;
1135:                 t[6][2] =  cx;
1136:                 t[4][3] =  cy;
1137:                 t[5][3] = -cx;
1138:              }
1139:              else {                               /* space frames */
1140:                 for(i = 4; i <= 6; i++)
1141:                     t[i][i] = 1.0;
1142:                 t[4][2] = -cz;
1143:                 t[4][3] =  cy;
1144:                 t[5][1] =  cz;
1145:                 t[5][3] = -cx;
1146:                 t[6][1] = -cy;
1147:                 t[6][2] =  cx;
1148:              }
1149:      }
1150: 
1151:      return(t);
1152: }
1153: 
1154: 
1155: /* 
1156:  *  =====================================================================
1157:  *  Modify Transformation_Matrix(T) : makes rows and columns to zero for
1158:  *  dof not connected to rigid body. Diagonal terms are retained as unity 
1159:  *  
1160:  *  Input - frame pointer frp, prob case  
1161:  *  =====================================================================
1162:  */ 
1163: 
1164: #ifdef __STDC__
1165: double **Modify_T_Matrix(double **T, EFRAME *frp, int no_dof, int rb_no)
1166: #else
1167: double **Modify_T_Matrix(T,frp, no_dof,rb_no)
1168: double **T;
1169: EFRAME *frp;
1170: int    no_dof, rb_no;
1171: #endif
1172: {
1173: int nn, i,j,k,ndof;
1174: int rb_elmt_type;
1175: 
1176:     rb_elmt_type  = frp->rigid[rb_no -1].rb_type;
1177: 
1178:     /* For Storey type rigid body                     */
1179:     /* all dofs(=3) assumed to be connected to R.Body */
1180: 
1181:     if(rb_elmt_type == 7) { /* Storey type rigid body ;dofs=3 */ 
1182:                             /* all dofs are connected hence returned unaltered     */
1183:                             /* printf("** In Modify T_Mat: Storey Type   \n");  */
1184:        return(T);
1185:     }
1186: 
1187:     /* general case */
1188: 
1189:     for(j = 1; j<= no_dof;j++){
1190:         ndof = frp->rigid[rb_no -1].rest_dof[j];
1191:         if(ndof >= 1) { /* dofs not connected by Rigid body are = 1 */
1192:            for(k = 1; k <= no_dof; k++)
1193:                T[j][k] =0.0;
1194:                T[j][j] = 1.0;
1195:         }
1196:     }
1197: 
1198:     return(T);
1199: }
1200: 
1201: /* 
1202:  *  ============================================
1203:  *  Set Element Properties for Feap problem        
1204:  *  Input   - p array                               
1205:  *  output  - p array                               
1206:  *  ============================================
1207:  */ 
1208: 
1209: #ifdef __STDC__
1210: ARRAY *Element_Property(ARRAY *p)
1211: #else
1212: ARRAY *Element_Property(p)
1213: ARRAY *p;
1214: #endif
1215: {
1216:    p = elmlib(p,PROPTY);
1217:    return(p);
1218: }
1219: 
1220: /* ================================================ */
1221: /* Set Material Properties for Feap problem         */
1222: /* Input  - p array                                 */
1223: /* output  - p array                                */
1224: /* ================================================ */
1225: 
1226: #ifdef __STDC__
1227: ARRAY *Mate_Property(ARRAY *p)
1228: #else
1229: ARRAY *Mate_Property(p)
1230: ARRAY *p;
1231: #endif
1232: {
1233: int isw = 15;
1234: 
1235:     p = elmlib(p,isw);
1236:     return(p);
1237: }
1238: 
1239: /* ================================================== */
1240: /* Assemble Element Property Vector                   */
1241: /* Input  - frame pointer frp, prob case              */
1242: /* ================================================== */
1243: 
1244: #ifdef __STDC__
1245: ARRAY *Eload_Property(ARRAY *p)
1246: #else
1247: ARRAY *Eload_Property(p)
1248: ARRAY *p;
1249: #endif
1250: {
1251: int  isw = 16;
1252: 
1253:     p = elmlib(p,isw) ;
1254:     return(p);
1255: }
1256:  
1257: #define PLASTIC 2
1258: 
1259: /* ======================================================*/
1260: /* function  double   *Transform_Force(frp,F,nn )        */
1261: /*  Calculates  Transform_Force Matrix  Ti   for         */
1262: /*  Relation:        Abi = Ti . Ami                      */
1263: /*  Converts the actions in Ami(at end of the member)    */
1264: /*  to the statically equivalent actions in Abi          */
1265: /*         ( at the working points of the rigid bodies ) */
1266: /* Input - frame pointer frp,                            */
1267: /*      F: Force Vector for the problem which is partially */
1268: /*            modified for node nu 'nn' in transforming  */
1269: /*            node force from this point to wkg pt of RB */
1270: /* ===================================================== */
1271: 
1272: #ifdef __STDC__
1273: QUANTITY *Transform_Force(EFRAME *frp, QUANTITY *F, int nn )
1274: #else
1275: QUANTITY *Transform_Force(frp,F,nn )
1276: EFRAME   *frp;
1277: QUANTITY *F;
1278: int nn ;
1279: #endif
1280: {
1281: int      i,k,rb_no,rb_elmt_type;
1282: double   **Fj, **Ff, **Tfj;
1283: double   cgx,cgy,cgz ,temp;
1284: 
1285:    k   = frp->no_dof;       /*for general case*/
1286:    Fj  = MatrixAllocIndirectDouble(k,1);
1287:    Tfj = MatrixAllocIndirectDouble(k,k);
1288: 
1289:    rb_no = frp->node[nn -1].rb_num;
1290:    rb_elmt_type = 7;
1291:    rb_elmt_type = frp->rigid[rb_no -1].rb_type;
1292: 
1293:    /* ------------------------------------------------------------------- */
1294:    /* compute coordinates of node wrt ref point- cg(wkg pt) of rigid body */
1295:    /* ------------------------------------------------------------------- */
1296: 
1297:     cgx = frp->node[nn-1].coord[1].value - frp->rigid[rb_no -1].xcg.value;
1298:     cgy = frp->node[nn-1].coord[2].value - frp->rigid[rb_no -1].ycg.value;
1299: 
1300: /********
1301:     if(frp->ndm >=3)
1302:        cgz  =   frp->node[nn-1].coord[3] - frp->rigid[rb_no -1].zcg;
1303:     else
1304:        cgz  = 0.0;
1305: 
1306: **********/
1307:     /* Transformation Matrix */
1308: 
1309:     Tfj = (double **)Transformation_Matrix(Tfj,(char *)rb_elmt_type,frp->no_dof, cgx,cgy,cgz);
1310:     Tfj = (double **)Modify_T_Matrix(Tfj, frp,frp->no_dof ,rb_no);
1311: 
1312:     /* Transformation Matrix- store for further use */
1313: 
1314:     frp->node[nn -1].TrT->uMatrix.daa = (double **) dMatrixTranspose(Tfj,frp->no_dof,frp->no_dof);
1315: 
1316:     for(i = 1; i<=k;i++)
1317:         Fj[i][1] = F[(nn  - 1) * frp->no_dof +i].value;
1318: 
1319:     Ff = (double **) dMatrixMult(Tfj, frp->no_dof, frp->no_dof,Fj, frp->no_dof, 1);
1320:     for(i = 1; i<=k;i++)
1321:         F[(nn  - 1) * frp->no_dof+i].value = Ff[i][1];
1322: 
1323:     MatrixFreeIndirectDouble(Ff, frp->no_dof);
1324:     MatrixFreeIndirectDouble(Fj, frp->no_dof);
1325:     MatrixFreeIndirectDouble(Tfj, frp->no_dof);
1326: 
1327:     return(F);
1328: }
1329: 
1330: 
1331: /* 
1332:  *  ============================================================================ 
1333:  *  Bound_Disp : Checks for displacement at restraint boundary nodes for an elmt             
1334:  *  
1335:  *  Input - frame pointer frp, element no  
1336:  *  ============================================================================ 
1337:  */ 
1338: 
1339: #ifdef __STDC__
1340: int Bound_Disp(EFRAME *frp, ARRAY *p, int elmt_no)
1341: #else
1342: int Bound_Disp(frp,p, elmt_no)
1343: EFRAME  *frp;
1344: ARRAY   *p;
1345: int     elmt_no;
1346: #endif
1347: {
1348: ELEMENT *el;
1349: double  displ;
1350: int     i,j,k,nen, no_dof,out;
1351: 
1352:     el  = &frp->element[elmt_no -1];      /* element ptr   */
1353:     p   = Assign_p_Array(frp,elmt_no,p, 10);
1354:     nen = p->nodes_per_elmt;
1355:     no_dof = p->dof_per_node;   /*ed feb 6  MIN(p->dof_per_node,frp->no_dof)*/
1356: 
1357:     out = NO;         /* default outcome  initialized */
1358: 
1359:     for(i = 1; i<=nen;i++){
1360:         for(j = 1; j<=no_dof;j++){
1361:             k = el->d_array[(i-1) *no_dof + j];
1362:             if(k<0) {           /* fixed or restrained dofs */
1363:                displ = frp->node[el->node_connect[i]  - 1].disp[j].value; 
1364:                p->nodal_loads[(i-1) * no_dof + j].value = displ;
1365:                if(displ != 0.0) out = YES;
1366:             }
1367:             else
1368:                p->nodal_loads[(i-1) * no_dof + j].value = 0.0;
1369:         }
1370:     }
1371: 
1372:     return(out);
1373: }
1374: 
1375: /* =========== */
1376: /* Modify_Load */
1377: /* =========== */
1378: 
1379: #ifdef __STDC__
1380: QUANTITY *Modify_Load(QUANTITY *b, EFRAME *frp, ARRAY *p, int elmt_no)
1381: #else
1382: QUANTITY *Modify_Load(b, frp, p, elmt_no)
1383: QUANTITY *b;
1384: EFRAME   *frp;
1385: ARRAY    *p;
1386: int      elmt_no;
1387: #endif
1388: {
1389: ELEMENT  *el;
1390: MATRIX   *Ke;
1391: int      *Ld;
1392: int      i,j,k;
1393: 
1394:     el = &frp->element[elmt_no -1];
1395:     p  = Assign_p_Array(frp,elmt_no,p, LOAD_MATRIX);
1396:     Ke = Element_Matrix(p,STIFF);
1397:     Ld = el->d_array;
1398: 
1399:     for(i = 1; i<=p->size_of_stiff;i++) {
1400:         k = Ld[i];
1401:         if(k>0) {
1402:            for(j = 1; j<=p->size_of_stiff;j++)
1403:                b[k].value = b[k].value - Ke->uMatrix.daa[i][j] * p->nodal_loads[j].value;
1404:         }
1405:     }
1406: 
1407:     return(b);
1408: }
1409: 
1410: 
1411: /* ============================================== */
1412: /*   Function to Add input loadvector pl          */
1413: /*                      to loadvector fvector     */
1414: /* ============================================== */
1415: 
1416: #ifdef __STDC__
1417: QUANTITY *Addload_Vector(EFRAME *frp, ARRAY *p, QUANTITY *fvector,  QUANTITY *pl, int elmt_no)
1418: #else
1419: QUANTITY *Addload_Vector(frp,p,fvector, pl, elmt_no)
1420: EFRAME   *frp;