1142
1143
1144
1146
1147
1148
1149#include "implicit_f.inc"
1150
1151
1152
1153#include "task_c.inc"
1154#include "param_c.inc"
1155#include "scr07_c.inc"
1156
1157
1158
1159 INTEGER INRBE3(*),ILRBE3(*),NG, NS,IROT,ID
1160
1162 . xyz(3,*), frbe3(6,*), skew(lskew,*),fdstnb(3,6,*), mdstnb(3,6,*)
1163
1164
1165
1166 INTEGER I, J, K,N, M ,NML, IAD,JJ,KG,NSNGLR,IELSUB,IERR,ng1
1167
1169 * tw(3,ng), rw(3,ng),
1170 * fufxlc(3,ng), fufylc(3,ng), fufzlc(3,ng),
1171 * fumxlc(3,ng), fumylc(3,ng), fumzlc(3,ng),
1172 * mxlc(3,ng), mylc(3,ng), mzlc(3,ng),
1173 * fufx(3,ng), fufy(3,ng), fufz(3,ng),
1174 * mufx(3,ng), mufy(3,ng), mufz(3,ng),
1175 * fumx(3,ng), fumy(3,ng), fumz(3,ng),
1176 * mx(3,ng), my(3,ng), mz(3,ng),
1177 * mumx(3,ng), mumy(3,ng), mumz(3,ng),
1178 * el(3,3,ng)
1180 * denfx, denfy, denfz, denmx, denmy, denmz,
1181 * refpt(3), cgmx(3), cgmy(3), cgmz(3), averef,
1182 * tfufx(3), tfufy(3), tfufz(3),
1183 * tmufx(3), tmufy(3), tmufz(3),
1184 * tfumx(3), tfumy(3), tfumz(3),
1185 * tmumx(3), tmumy(3), tmumz(3),
1186 * a(6,6), c(6,6), t(3,3)
1187
1188
1189
1190 IF (ng==0) RETURN
1191 CALL zero1(fdstnb,3*ng*6)
1192 IF (irot>0)
CALL zero1(mdstnb,3*ng*6)
1198 ierr = 0
1199
1200 refpt(1) = xyz(1,ns)
1201 refpt(2) = xyz(2,ns)
1202 refpt(3) = xyz(3,ns)
1203 DO k = 1, ng
1204 DO i = 1, 3
1205 tw(i,k) = frbe3(i,k)
1206 rw(i,k) = frbe3(i+3,k)
1207 ENDDO
1208 ENDDO
1209
1210
1211
1212
1213
1214 IF (ng == 2.AND.irot==0) THEN
1215 ierr = 322
1216 GOTO 999
1217 ENDIF
1218
1219
1220
1221 DO k = 1, ng
1222 ielsub = ilrbe3(k)
1223 IF (ielsub > 0) THEN
1224 DO i = 1, 3
1225 el(i,1,k) = skew(i,ielsub)
1226 el(i,2,k) = skew(i+3,ielsub)
1227 el(i,3,k) = skew(i+6,ielsub)
1228 ENDDO
1229 ENDIF
1230 ENDDO
1231
1232
1233
1234 denfx = zero
1235 denfy = zero
1236 denfz = zero
1237 averef = zero
1238
1239 DO 70 k = 1, ng
1240 kg = inrbe3(k)
1241 ielsub = ilrbe3(k)
1242 IF (ielsub > 0) THEN
1243
1244
1245
1246 DO 60 i = 1, 3
1247 denfx = denfx + tw(i,k)*el(i,1,k)**2
1248 denfy = denfy + tw(i,k)*el(i,2,k)**2
1249 denfz = denfz + tw(i,k)*el(i,3,k)**2
1250 60 CONTINUE
1251 ELSE
1252 denfx = denfx + tw(1,k)
1253 denfy = denfy + tw(2,k)
1254 denfz = denfz + tw(3,k)
1255 END IF
1256
1257 averef = averef + sqrt( (xyz(1,kg) - refpt(1))**2 +
1258 * (xyz(2,kg) - refpt(2))**2 +
1259 * (xyz(3,kg) - refpt(3))**2 )
1260 70 CONTINUE
1261
1262 IF (abs(denfx) <= em20) THEN
1263 ierr = 326
1264 ENDIF
1265
1266 IF (abs(denfy) <= em20) THEN
1267 ierr = 327
1268 ENDIF
1269
1270 IF (abs(denfz) <= em20) THEN
1271 ierr = 328
1272 ENDIF
1273 IF (ierr /= 0) GOTO 999
1274 averef = averef/ng
1275 IF (averef == zero) averef = one
1276
1277
1278
1279
1280 DO 40 k = 1, ng
1281 kg = inrbe3(k)
1282 ielsub = ilrbe3(k)
1283 IF (ielsub > 0) THEN
1284
1285
1286
1287 DO 10 i = 1, 3
1288 cgmx(2) = cgmx(2) + tw(i,k)*el(i,3,k)**2*xyz(2,kg)
1289 cgmx(3) = cgmx(3) + tw(i,k)*el(i,2,k)**2*xyz(3,kg)
1290 10 CONTINUE
1291
1292 DO 20 i = 1, 3
1293 cgmy(3) = cgmy(3) + tw(i,k)*el(i,1,k)**2*xyz(3,kg)
1294 cgmy(1) = cgmy(1) + tw(i,k)*el(i,3,k)**2*xyz(1,kg)
1295 20 CONTINUE
1296
1297 DO 30 i = 1, 3
1298 cgmz(1) = cgmz(1) + tw(i,k)*el(i,2,k)**2*xyz(1,kg)
1299 cgmz(2) = cgmz(2) + tw(i,k)*el(i,1,k)**2*xyz(2,kg)
1300 30 CONTINUE
1301
1302 ELSE
1303 cgmx(2) = cgmx(2) + tw(3,k)*xyz(2,kg)
1304 cgmx(3) = cgmx(3) + tw(2,k)*xyz(3,kg)
1305
1306 cgmy(3) = cgmy(3) + tw(1,k)*xyz(3,kg)
1307 cgmy(1) = cgmy(1) + tw(3,k)*xyz(1,kg)
1308
1309 cgmz(1) = cgmz(1) + tw(2,k)*xyz(1,kg)
1310 cgmz(2) = cgmz(2) + tw(1,k)*xyz(2,kg)
1311 END IF
1312 40 CONTINUE
1313 cgmx(2) = cgmx(2)/denfz
1314 cgmx(3) = cgmx(3)/denfy
1315
1316 cgmy(3) = cgmy(3)/denfx
1317 cgmy(1) = cgmy(1)/denfz
1318
1319 cgmz(1) = cgmz(1)/denfy
1320 cgmz(2) = cgmz(2)/denfx
1321
1322 denmx = zero
1323 denmy = zero
1324 denmz = zero
1325
1326 DO 90 k = 1, ng
1327 kg = inrbe3(k)
1328 ielsub = ilrbe3(k)
1329
1330
1331
1332
1333
1334
1335 IF (ielsub > 0) THEN
1336
1337
1338
1339 DO 80 i = 1, 3
1340 denmx = denmx + rw(i,k)*el(i,1,k)**2*averef**2 +
1341 * tw(i,k)*( el(i,3,k)*(xyz(2,kg) - cgmx(2)) -
1342 * el(i,2,k)*(xyz(3,kg) - cgmx(3))
1343 * ) **2
1344 denmy = denmy + rw(i,k)*el(i,2,k)**2*averef**2 +
1345 * tw(i,k)*( el(i,1,k)*(xyz(3,kg) - cgmy(3)) -
1346 * el(i,3,k)*(xyz(1,kg) - cgmy(1))
1347 * ) **2
1348 denmz = denmz + rw(i,k)*el(i,3,k)**2*averef**2 +
1349 * tw(i,k)*( el(i,2,k)*(xyz(1,kg) - cgmz(1)) -
1350 * el(i,1,k)*(xyz(2,kg) - cgmz(2))
1351 * ) **2
1352 80 CONTINUE
1353 ELSE
1354 denmx = denmx + rw(1,k)*averef**2 +
1355 * tw(2,k)*(xyz(3,kg) - cgmx(3))**2 +
1356 * tw(3,k)*(xyz(2,kg) - cgmx(2))**2
1357 denmy = denmy + rw(2,k)*averef**2 +
1358 * tw(1,k)*(xyz(3,kg) - cgmy(3))**2 +
1359 * tw(3,k)*(xyz(1,kg) - cgmy(1))**2
1360 denmz = denmz + rw(3,k)*averef**2 +
1361 * tw(2,k)*(xyz(1,kg) - cgmz(1))**2 +
1362 * tw(1,k)*(xyz(2,kg) - cgmz(2))**2
1363 END IF
1364 90 CONTINUE
1365
1366
1367
1368
1369
1370 IF (abs(denmx) <= em20) THEN
1371 ierr = 329
1372 ENDIF
1373
1374 IF (abs(denmy) <= em20) THEN
1375 ierr = 330
1376 ENDIF
1377
1378 IF (abs(denmz) <= em20) THEN
1379 ierr = 331
1380 ENDIF
1381
1382 IF (ierr /= 0) GOTO 999
1383
1384
1385
1386
1387 CALL rbe3uf(inrbe3,ilrbe3,el,tw,xyz,refpt,
1388 * fufxlc,fufylc,fufzlc,fufx,fufy,fufz,mufx,mufy,mufz,
1389 * tfufx,tfufy,tfufz,tmufx,tmufy,tmufz,
1390 * denfx,denfy,denfz,ng)
1391
1392
1393
1394
1395
1396 CALL rbe3um(inrbe3,ilrbe3,el,tw,rw,xyz,refpt,cgmx,cgmy,cgmz,
1397 * fumxlc,fumylc,fumzlc,mxlc,mylc,mzlc,
1398 * fumx,fumy,fumz,mx,my,mz,mumx,mumy,mumz,
1399 * tfumx,tfumy,tfumz,tmumx,tmumy,tmumz,
1400 * averef,denmx,denmy,denmz,ng,irot )
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422 DO 120 i = 1, 3
1423 k = i + 3
1424 a(i,1) = tfufx(i)
1425 a(k,1) = tmufx(i)
1426 a(i,2) = tfufy(i)
1427 a(k,2) = tmufy(i)
1428 a(i,3) = tfufz(i)
1429 a(k,3) = tmufz(i)
1430 a(i,4) = tfumx(i)
1431 a(k,4) = tmumx(i)
1432 a(i,5) = tfumy(i)
1433 a(k,5) = tmumy(i)
1434 a(i,6) = tfumz(i)
1435 a(k,6) = tmumz(i)
1436 120 CONTINUE
1437
1438
1439
1440 nsnglr = 0
1441 CALL invert(a,c,6,nsnglr)
1442 IF (nsnglr /= 0) THEN
1443 ierr = 332
1444 GOTO 999
1445 ENDIF
1446
1447 DO i = 1, 3
1448 DO j = 1, 6
1449 DO k = 1, ng
1450 fdstnb(i,j,k) = c(1,j)*fufx(i,k) + c(2,j)*fufy(i,k) +
1451 * c(3,j)*fufz(i,k) + c(4,j)*fumx(i,k) +
1452 * c(5,j)*fumy(i,k) + c(6,j)*fumz(i,k)
1453 ENDDO
1454 ENDDO
1455 ENDDO
1456 IF (irot>0) THEN
1457 DO i = 1, 3
1458 DO j = 1, 6
1459 DO k = 1, ng
1460 mdstnb(i,j,k) = c(4,j)*mx(i,k) + c(5,j)*my(i,k) +
1461 * c(6,j)*mz(i,k)
1462 ENDDO
1463 ENDDO
1464 ENDDO
1465 END IF
1466
1467 999 CONTINUE
1468 IF (ierr>0) THEN
1469 IF(ispmd==0)THEN
1470 CALL ancmsg(msgid=108,anmode=aninfo,
1472 ENDIF
1473 mstop = 1
1474 ENDIF
1475
1476
1477
1478 RETURN
subroutine invert(matrix, inverse, n, errorflag)
subroutine rbe3um(inrbe3, ilrbe3, el, tw, rw, xyz, refpt, cgmx, cgmy, cgmz, fumxlc, fumylc, fumzlc, mxlc, mylc, mzlc, fumx, fumy, fumz, mx, my, mz, mumx, mumy, mumz, tfumx, tfumy, tfumz, tmumx, tmumy, tmumz, averef, denmx, denmy, denmz, ng, irot)
subroutine rbe3uf(inrbe3, ilrbe3, el, tw, xyz, refpt, fufxlc, fufylc, fufzlc, fufx, fufy, fufz, mufx, mufy, mufz, tfufx, tfufy, tfufz, tmufx, tmufy, tmufz, denfx, denfy, denfz, ng)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)