428
429
430
433
434
435
436#include "implicit_f.inc"
437#include "comlock.inc"
438
439
440
441#include "mvsiz_p.inc"
442
443
444
445#include "scr03_c.inc"
446#include "scr07_c.inc"
447#include "scr14_c.inc"
448#include "scr16_c.inc"
449#include "com06_c.inc"
450#include "com08_c.inc"
451#include "parit_c.inc"
452#include "scr18_c.inc"
453#include "impl1_c.inc"
454
455
456
457 INTEGER (*), NTY, IMAST, LFT, LLT, NFT
459 INTEGER IRECT(4,*), MSR(*), NSV(*), IRTL(*), IRTLO(*), ISKY(*)
460 my_real x(3,*), e(*), cst(2,*), fric0(3,*), fsav(*),
461 . fskyi(lskyi,nfskyi),
462 . v(3,*), fcont(3,*),cf(*), freq, frot_p(*),ftsav(*),
463 . ftcont(3,*)
464 TYPE(H3D_DATABASE) ::
465 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX1,IX2,IX3,IX4
466 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: n1,n2,n3
467 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: xp,yp,zp
468 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: ssc,ttc,xface,stif
469 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: h1,h2,h3,h4,
area,fni
470
471
472
473 INTEGER I, IL, LOLD, JJ, NN, J3,J2, J1, IG, I3, I2, I1, K, IFQ, MFROT, NISKYL
475 . h(4), xx1(4), xx2(4), xx3(4),
476 . fxi(mvsiz), fyi(mvsiz),
477 . fzi(mvsiz), fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz), fy1(mvsiz), fy2(mvsiz),
478 . fy3(mvsiz), fy4(mvsiz), fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
479 . ss0, tt0, xc, econvt,
alpha, alphi,
480 . yc, zc, xc0, yc0, zc0, sp, sm, tp, tm, ansx, ansy, ansz, fmax,
481 . stf, fti, fn, tn1, tn2, tn3, tn, dtm, xmu, vx,vy,vz,vv,v2,p,
482 . vv1,vv2,v21,dmu,aa
483 INTEGER IRTLO_I(MVSIZ)
485
486
487
488
489 IF (inconv/=1) THEN
490 DO i=lft,llt
491 il=i+nft
492 irtlo_i(i)=irtlo(il)
493 fric0_i(1,i)=fric0(1,il)
494 fric0_i(2,i)=fric0(2,il)
495 fric0_i(3,i)=fric0(3,il)
496 ENDDO
497 ENDIF
498 DO 300 i=lft,llt
499 fxi(i)=zero
500 fyi(i)=zero
501 fzi(i)=zero
502 il=i+nft
503 econvt = zero
504 IF(xface(i)==zero) THEN
505
506
507
508 irtlo(il)=0
509 fric0(1,il)=zero
510 fric0(2,il)=zero
511 fric0(3,il)=zero
512 ELSE
513
514 lold=iabs(irtlo(il))
515 IF(lold==0)THEN
516
517
518
519 irtlo(il)=irtl(il)*xface(i)
520 cst(1,il)=ssc(i)
521 cst(2,il)=ttc(i)
522 ELSE
523
524
525
526 ss0=cst(1,il)
527 tt0=cst(2,il)
528 fxi(i)=fric0(1,il)
529 fyi(i)=fric0(2,il)
530 fzi(i)=fric0(3,il)
531
532 xc=xp(i)
533 yc=yp(i)
534 zc=zp(i)
535 DO jj=1,4
536 nn=msr(irect(jj,lold))
537 xx1(jj)=x(1,nn)
538 xx2(jj)=x(2,nn)
539 xx3(jj)=x(3,nn)
540 ENDDO
541 xc0=zero
542 yc0=zero
543 zc0=zero
544 sp=one+ss0
545 sm=one-ss0
546 tp=fourth*(one + tt0)
547 tm=fourth*(one - tt0)
548 h(1)=tm*sm
549 h(2)=tm*sp
550 h(3)=tp*sp
551 h(4)=tp*sm
552 DO jj=1,4
553 xc0=xc0+h(jj)*xx1(jj)
554 yc0=yc0+h(jj)*xx2(jj)
555 zc0=zc0+h(jj)*xx3(jj)
556 ENDDO
557 ansx= (xc-xc0)
558 ansy= (yc-yc0)
559 ansz= (zc-zc0)
560 econvt = econvt + half*(fxi(i)*ansx+fyi(i)*ansy+fzi(i)*ansz)
561 stf=em01*stif(i)
562 fxi(i)=fxi(i) + ansx*stf
563 fyi(i)=fyi(i) + ansy*stf
564 fzi(i)=fzi(i) + ansz*stf
565
566 IF (codvers>=44) THEN
567 ifq = ipari(31)
568 IF (ifq>0) THEN
569 IF (ifq==3) freq =
max(one,freq*dt12)
572 k = 3*(il-1)
573 IF (fni(i)==0) THEN
574 ftsav(k+1) = zero
575 ftsav(k+2) = zero
576 ftsav(k+3) = zero
577 ELSE
578 fxi(i)=
alpha*fxi(i) + alphi*ftsav(k+1)
579 fyi(i)=
alpha*fyi(i) + alphi*ftsav(k+2)
580 fzi(i)=
alpha*fzi(i) + alphi*ftsav(k+3)
581 IF (inconv==1) THEN
582 ftsav(k+1) = fxi(i)
583 ftsav(k+2) = fyi(i)
584 ftsav(k+3) = fzi(i)
585 ENDIF
586 ENDIF
587 ENDIF
588 ENDIF
589 fti=sqrt(fxi(i)*fxi(i)+fyi(i)*fyi(i)+fzi(i)*fzi(i))
590
591 fn =fxi(i)*n1(i)+fyi(i)*n2(i)+fzi(i)*n3(i)
592 tn1=fxi(i)-n1(i)*fn
593 tn2=fyi(i)-n2(i)*fn
594 tn3=fzi(i)-n3(i)*fn
595 tn =sqrt(tn1*tn1+tn2*tn2+tn3*tn3)
596 IF(tn/=zero)THEN
597 tn1=tn1/tn
598 tn2=tn2/tn
599 tn3=tn3/tn
600 ELSE
601 tn3=zero
602 tn=sqrt(n1(i)*n1(i)+n2(i)*n2(i))
603 IF(tn/=zero)THEN
604 tn2=-n1(i)/tn
605 tn1= n2(i)/tn
606 ELSE
607 tn2=zero
608 tn1=one
609 ENDIF
610 ENDIF
611
613 xmu = fric
614 IF(codvers>=44) THEN
615
616
617
618 mfrot = ipari(30)
619 IF(mfrot>=1) THEN
620 ig=nsv(il)
621 vx = v(1,ig) - h(1)*v(1,ix1(i)) - h(2)*v(1,ix2(i))
622 . - h(3)*v(1,ix3(i)) - h(4)*v(1,ix4(i))
623 vy = v(2,ig) - h(1)*v(2,ix1(i)) - h(2)*v(2,ix2(i))
624 . - h(3)*v(2,ix3(i)) - h(4)*v(2,ix4(i))
625 vz = v(3,ig) - h(1)*v(3,ix1(i)) - h(2)*v(3,ix2(i))
626 . - h(3)*v(3,ix3(i)) - h(4)*v(3,ix4(i))
627
628 aa = n1(i)*vx + n2(i)*vy + n3(i)*vz
629 vx = vx - n1(i)*aa
630 vy = vy - n2(i)*aa
631 vz = vz - n3(i)*aa
632 v2 = vx*vx + vy*vy + vz*vz
633 vv = sqrt(
max(em30, v2))
634 ENDIF
635 IF(mfrot==1)THEN
636 xmu = fric + (frot_p(1) + frot_p(4)*p ) * p
637 . +(frot_p(2) + frot_p(3)*p) * vv + frot_p(5)*v2
638 ELSEIF(mfrot==2)THEN
639
640 xmu = fric
641 . + frot_p(1)*exp(frot_p(2)*vv)*p**2
642 . + frot_p(3)*exp(frot_p(4)*vv)*p
643 . + frot_p(5)*exp(frot_p(6)*vv)
644 ELSEIF(mfrot==3)THEN
645
646 IF(vv>=0.AND.vv<=frot_p(5)) THEN
647 dmu = frot_p(3)-frot_p(1)
648 vv1 = vv / frot_p(5)
649 xmu = frot_p(1)+ dmu*vv1*(two-vv1)
650 ELSEIF(vv>frot_p(5).AND.vv<frot_p(6)) THEN
651 dmu = frot_p(4)-frot_p(3)
652 vv1 = (vv - frot_p(5))/(frot_p(6)-frot_p(5))
653 xmu
654 ELSE
655 dmu = frot_p(2)-frot_p(4)
656 vv2 = (vv - frot_p(6))**2
657 xmu = frot_p(2) - dmu / (one + dmu*vv2)
658 ENDIF
659 ENDIF
660 fmax = -
min(xmu*fni(i),zero)
661 IF(fti>fmax)THEN
662
663 fxi(i)=tn1*fmax
664 fyi(i)=tn2*fmax
665 fzi(i)=tn3*fmax
666 fric0(1,il)=fxi(i)
667 fric0(2,il)=fyi(i)
668 fric0(3,il)=fzi(i)
669 irtlo(il)=irtl(il)*xface(i)
670 cst(1,il)=ssc(i)
671 cst(2,il)=ttc(i)
672 ELSE
673
674 fxi(i)=tn1*fti
675 fyi(i)=tn2*fti
676 fzi(i)=tn3*fti
677 ENDIF
678 ELSE
679
680
681
682 ig=nsv(il)
683 vx = v(1,ig) - h(1)*v(1,ix1(i)) - h(2)*v(1,ix2(i))
684 . - h(3)*v(1,ix3(i)) - h(4)*v(1,ix4(i))
685 vy = v(2,ig) - h(1)*v(2,ix1(i)) - h(2)*v(2,ix2(i))
686 . - h(3)*v(2,ix3(i)) - h(4)*v(2,ix4(i))
687 vz = v(3,ig) - h(1)*v(3,ix1(i)) - h(2)*v(3,ix2(i))
688 . - h(3)*v(3,ix3(i)) - h(4)*v(3,ix4(i))
689 vv = vx*vx + vy*vy + vz*vz
690 xmu = fric + ( cf(1) + cf(4)*p ) * p
691 . + ( cf(2) + cf(3)*p ) * sqrt(vv) + cf(5)*vv
692 fmax= -
min(xmu*fni(i),zero)
693 IF(fti>fmax)THEN
694
695
696
697 fxi(i)=tn1*fmax
698 fyi(i)=tn2*fmax
699 fzi(i)=tn3*fmax
700 fric0(1,il)=fxi(i)
701 fric0(2,il)=fyi(i)
702 fric0(3,il)=fzi(i)
703 irtlo(il)=irtl(il)*xface(i)
704 cst(1,il)=ssc(i)
705 cst(2,il)=ttc(i)
706 ELSE
707
708
709
710 fxi(i)=tn1*fti
711 fyi(i)=tn2*fti
712 fzi(i)=tn3*fti
713 ENDIF
714 ENDIF
715 econvt = econvt + half*(fxi(i)*ansx+fyi(i)*ansy+fzi(i)*ansz)
716 ENDIF
717 ENDIF
718
719 300 CONTINUE
720
721 IF (inconv==1) THEN
722
723 econtv = econtv + econvt
724
725 fsav(27) = fsav(27) + econvt
726
727
728
729 dtm=imast*dt12
730 DO i=lft,llt
731 fsav(4)=fsav(4)+fxi(i)*dtm
732 fsav(5)=fsav(5)+fyi(i)*dtm
733 fsav(6)=fsav(6)+fzi(i)*dtm
734 ENDDO
735 END IF
736
737 DO i=lft,llt
738 fx1(i)=fxi(i)*h1(i)
739 fy1(i)=fyi(i)*h1(i)
740 fz1(i)=fzi(i)*h1(i)
741
742 fx2(i)=fxi(i)*h2(i)
743 fy2(i)=fyi(i)*h2(i)
744 fz2(i)=fzi(i)*h2(i)
745
746 fx3(i)=fxi(i)*h3(i)
747 fy3(i)=fyi(i)*h3(i)
748 fz3(i)=fzi(i)*h3(i)
749
750 fx4(i)=fxi(i)*h4(i)
751 fy4(i)=fyi(i)*h4(i)
752 fz4(i)=fzi(i)*h4(i)
753 ENDDO
754
755 IF(iparit==0)THEN
756
757
758
759 DO i=lft,llt
760 il=i+nft
761 j3=3*ix1(i)
762 j2=j3-1
763 j1=j2-1
764 e(j1)=e(j1)+fx1(i)
765 e(j2)=e(j2)+fy1(i)
766 e(j3)=e(j3)+fz1(i)
767
768 j3=3*ix2(i)
769 j2=j3-1
770 j1=j2-1
771 e(j1)=e(j1)+fx2(i)
772 e(j2)=e(j2)+fy2(i)
773 e(j3)=e(j3)+fz2(i)
774
775 j3=3*ix3(i)
776 j2=j3-1
777 j1=j2-1
778 e(j1)=e(j1)+fx3(i)
779 e(j2)=e(j2)+fy3(i)
780 e(j3)=e(j3)+fz3(i)
781
782 j3=3*ix4(i)
783 j2=j3-1
784 j1=j2-1
785 e(j1)=e(j1)+fx4(i)
786 e(j2)=e(j2)+fy4(i)
787 e(j3)=e(j3)+fz4(i)
788
789
790
791 ig=nsv(il)
792 i3=3*ig
793 i2=i3-1
794 i1=i2-1
795 e(i1)=e(i1)-fxi(i)
796 e(i2)=e(i2)-fyi(i)
797 e(i3)=e(i3)-fzi(i)
798
799 ENDDO
800
801 ELSE
802
803#include "lockon.inc"
804 niskyl = nisky
805 nisky = nisky + 5 * llt
806#include "lockoff.inc"
807
808 IF(kdtint==0)THEN
809 DO 440 i=lft,llt
810 niskyl = niskyl + 1
811 fskyi(niskyl,1)=fx1(i)
812 fskyi(niskyl,2)=fy1(i)
813 fskyi(niskyl,3)=fz1(i)
814 fskyi(niskyl,4)=zero
815 isky(niskyl) = ix1(i)
816 niskyl = niskyl + 1
817 fskyi(niskyl,1)=fx2(i)
818 fskyi(niskyl,2)=fy2(i)
819 fskyi(niskyl,3)=fz2(i)
820 fskyi(niskyl,4)=zero
821 isky(niskyl) = ix2(i)
822 niskyl = niskyl + 1
823 fskyi(niskyl,1)=fx3(i)
824 fskyi(niskyl,2)=fy3(i)
825 fskyi(niskyl,3)=fz3(i)
826 fskyi(niskyl,4)=zero
827 isky(niskyl) = ix3(i)
828 niskyl = niskyl + 1
829 fskyi(niskyl,1)=fx4(i)
830 fskyi(niskyl,2)=fy4(i)
831 fskyi(niskyl,3)=fz4(i)
832 fskyi(niskyl,4)=zero
833 isky(niskyl) = ix4(i)
834 niskyl = niskyl + 1
835 fskyi(niskyl,1)=-fxi(i)
836 fskyi(niskyl,2)=-fyi(i)
837 fskyi(niskyl,3)=-fzi(i)
838 fskyi(niskyl,4)=zero
839 il=i+nft
840 isky(niskyl) = nsv(il)
841 440 CONTINUE
842 ELSE
843 DO i=lft,llt
844 niskyl = niskyl + 1
845 fskyi(niskyl,1)=fx1(i)
846 fskyi(niskyl,2)=fy1(i)
847 fskyi(niskyl,3)=fz1(i)
848 fskyi(niskyl,4)=zero
849 fskyi(niskyl,5)=zero
850 isky(niskyl) = ix1(i)
851 niskyl = niskyl + 1
852 fskyi(niskyl,1)=fx2(i)
853 fskyi(niskyl,2)=fy2(i)
854 fskyi(niskyl,3)=fz2(i)
855 fskyi(niskyl,4)=zero
856 fskyi(niskyl,5)=zero
857 isky(niskyl) = ix2(i)
858 niskyl = niskyl + 1
859 fskyi(niskyl,1)=fx3(i)
860 fskyi(niskyl,2)=fy3(i)
861 fskyi(niskyl,3)=fz3(i)
862 fskyi(niskyl,4)=zero
863 fskyi(niskyl,5)=zero
864 isky(niskyl) = ix3(i)
865 niskyl = niskyl + 1
866 fskyi(niskyl,1)=fx4(i)
867 fskyi(niskyl,2)=fy4(i)
868 fskyi(niskyl,3)=fz4(i)
869 fskyi(niskyl,4)=zero
870 fskyi(niskyl,5)=zero
871 isky(niskyl) = ix4(i)
872 niskyl = niskyl + 1
873 fskyi(niskyl,1)=-fxi(i)
874 fskyi(niskyl,2)=-fyi(i)
875 fskyi(niskyl,3)=-fzi(i)
876 fskyi(niskyl,4)=zero
877 fskyi(niskyl,5)=zero
878 il=i+nft
879 isky(niskyl) = nsv(il)
880 ENDDO
881 ENDIF
882 ENDIF
883
884 IF (inconv/=1) THEN
885 DO i=lft,llt
886 il=i+nft
887 irtlo(il)=irtlo_i(i)
888 fric0(1,il)=fric0_i(1,i)
889 fric0(2,il)=fric0_i(2,i)
890 fric0(3,il)=fric0_i(3,i)
891 ENDDO
892 RETURN
893 ENDIF
894
895 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0 .AND.
896 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
897 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D /= 0))THEN
898#include "lockon.inc"
899 DO i=1,llt
900 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
901 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
902 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
903 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
904 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
905 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
906 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
907 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
908 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
909 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
910 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
911 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
912 fcont(1,nsv(i+nft))=fcont(1,nsv(i+nft))- fxi(i)
913 fcont(2,nsv(i+nft))=fcont(2,nsv(i+nft))- fyi(i)
914 fcont(3,nsv(i+nft))=fcont(3,nsv(i+nft))- fzi(i)
915 ENDDO
916#include "lockoff.inc"
917 ENDIF
918
919 IF(anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
920 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP) .OR.
921 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))THEN
922#include "lockon.inc"
923 DO i=1,llt
924 ftcont(1,ix1(i)) =ftcont(1,ix1(i)) + fx1(i)
925 ftcont(2,ix1(i)) =ftcont(2,ix1(i)) + fy1(i)
926 ftcont(3,ix1(i)) =ftcont(3,ix1(i)) + fz1(i)
927 ftcont(1,ix2(i)) =ftcont(1,ix2(i)) + fx2(i)
928 ftcont(2,ix2(i)) =ftcont(2,ix2(i)) + fy2(i)
929 ftcont(3,ix2(i)) =ftcont(3,ix2(i)) + fz2(i)
930 ftcont(1,ix3(i)) =ftcont(1,ix3(i)) + fx3(i)
931 ftcont(2,ix3(i)) =ftcont(2,ix3(i)) + fy3(i)
932 ftcont(3,ix3(i)) =ftcont(3,ix3(i)) + fz3(i)
933 ftcont(1,ix4(i)) =ftcont(1,ix4(i)) + fx4(i)
934 ftcont(2,ix4(i)) =ftcont(2,ix4(i)) + fy4(i)
935 ftcont(3,ix4(i)) =ftcont(3,ix4(i)) + fz4(i)
936 ftcont(1,nsv(i+nft))=ftcont(1,nsv(i+nft))- fxi(i)
937 ftcont(2,nsv(i+nft))=ftcont(2,nsv(i+nft))- fyi(i)
938 ftcont(3,nsv(i+nft))=ftcont(3,nsv(i+nft))- fzi(i)
939 ENDDO
940#include "lockoff.inc"
941 ENDIF
942
943 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)