403
404
405
406#include "implicit_f.inc"
407
408
409
410#include "param_c.inc"
411#include "com01_c.inc"
412
413
414
415
417 . tt,dt1,stif,vx,vy,vz,b1,b2,b3,c1,c2,c3,det
419 . skew(9),rs(3),rm(3),rx(4),ry(4),rz(4),
420 . va(3), vb(3), vc(3), vd(3),vrs(3),vrm(3)
421
422
423
424
426 . r(3),dwdu,
427 . x12,x22,x32,x42,y12,y22,y32,y42,z12,z22,z32,z42,
428 . xx,yy,zz,xxx,yyy,zzz,xy,yz,zx,xy2,yz2,zx2,
429 . bb1,bb2,bb3,cc1,cc2,cc3,
430 . mgx,mgy,mgz,wx,wy,wz,in_secnd
431
432
433 IF (tt == zero) THEN
434
435 skew(1) = rs(1)-rm(1)
436 skew(2) = rs(2)-rm(2)
437 skew(3) = rs(3)-rm(3)
438 ENDIF
439
440
441 x12=rx(1)*rx(1)
442 x22=rx(2)*rx(2)
443 x32=rx(3)*rx(3)
444 x42=rx(4)*rx(4)
445 y12=ry(1)*ry(1)
446 y22=ry(2)*ry(2)
447 y32=ry(3)*ry(3)
448 y42=ry(4)*ry(4)
449 z12=rz(1)*rz(1)
450 z22=rz(2)*rz(2)
451 z32=rz(3)*rz(3)
452 z42=rz(4)*rz(4)
453 xx=x12 + x22 + x32 + x42
454 yy=y12 + y22 + y32 + y42
455 zz=z12 + z22 + z32 + z42
456 xy=rx(1)*ry(1) + rx(2)*ry(2) + rx(3)*ry(3) + rx(4)*ry(4)
457 yz=ry(1)*rz(1) + ry(2)*rz(2) + ry(3)*rz(3) + ry(4)*rz(4)
458 zx=rz(1)*rx(1) + rz(2)*rx(2) + rz(3)*rx(3) + rz(4)*rx(4)
459 zzz=xx+yy
460 xxx=yy+zz
461 yyy=zz+xx
462 xy2=xy*xy
463 yz2=yz*yz
464 zx2=zx*zx
465 det= xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2
466 . - two*xy*yz*zx
467 det=one/det
468 b1=zzz*yyy-yz2
469 b2=xxx*zzz-zx2
470 b3=yyy*xxx-xy2
471 c3=zzz*xy+yz*zx
472 c1=xxx*yz+zx*xy
473 c2=yyy*zx+xy*yz
474
475 mgx = ry(1)*va(3) + ry(2)*vb(3) + ry(3)*vc(3) + ry(4)*vd(3)
476 . - rz(1)*va(2) - rz(2)*vb(2) - rz(3)*vc(2) - rz(4)*vd(2)
477 mgy = rz(1)*va(1) + rz(2)*vb(1) + rz(3)*vc(1) + rz(4)*vd(1)
478 . - rx(1)*va(3) - rx(2)*vb(3) - rx(3)*vc(3) - rx(4)*vd(3)
479 mgz = rx(1)*va(2) + rx(2)*vb(2) + rx(3)*vc(2) + rx(4)*vd(2)
480 . - ry(1)*va(1) - ry(2)*vb(1)
481
482 vrm(1)=det*(mgx*b1+mgy*c3+mgz*c2)
483 vrm(2)=det*(mgy*b2+mgz*c1+mgx*c3)
484 vrm(3)=det*(mgz*b3+mgx*c2+mgy*c1)
485
486 IF (iroddl == 0 .OR. in_secnd == zero) THEN
487
488 wx = vrm(1)
489 wy = vrm(2)
490 wz = vrm(3)
491 ELSE
492 wx = (vrm(1)+vrs(1))*half
493 wy = (vrm(2)+vrs(2))*half
494 wz = (vrm(3)+vrs(3))*half
495 ENDIF
496
497
498 r(1)=skew(1)
499 r(2)=skew(2)
500 r(3)=skew(3)
501
502
503 vx = vx - (wy*r(3)-wz*r(2)) -(vrm(2)*rm(3)-vrm(3)*rm(2))
504 vy = vy - (wz*r(1)-wx*r(3)) -(vrm(3)*rm(1)-vrm(1)*rm(3))
505 vz = vz - (wx*r(2)-wy*r(1)) -(vrm(1)*rm(2)-vrm(2)*rm(1))
506
507 bb1=b1*b1
508 bb2=b2*b2
509 bb3=b3*b3
510 cc1=c1*c1
511 cc2=c2*c2
512 cc3=c3*c3
513 dwdu=det*sqrt(
max(bb1*(yy+zz)+cc3*(zz+xx)+cc2*(xx+yy),
514 . bb2*(zz+xx)+cc1*(xx+yy)+cc3*(yy+zz),
515 . bb3*(xx+yy)+cc2*(yy+zz)+cc1*(zz+xx)))
516
517 stif=sqrt((r(1)*r(1)+r(2)*r(2)+r(3)*r(3)))*dwdu
518
519
520 RETURN