428
429
430
432 USE output_mod, ONLY: output_
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 TYPE(OUTPUT_), INTENT(inout) :: OUTPUT
458 INTEGER IPARI(*), NTY, IMAST, LFT, LLT, NFT
460 INTEGER IRECT(4,*), MSR(*), NSV(*), IRTL(*), IRTLO(*), ISKY(*)
461 my_real x(3,*), e(*), cst(2,*), fric0(3,*), fsav(*),
462 . fskyi(lskyi,nfskyi),
463 . v(3,*), fcont(3,*),cf(*), freq, frot_p(*),ftsav(*),
464 . ftcont(3,*)
465 TYPE(H3D_DATABASE) :: H3D_DATA
466 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX1,IX2,IX3,IX4
467 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: n1,n2,n3
468 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: xp,yp,zp
469 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: ssc,ttc,xface,stif
470 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: h1,h2,h3,h4,
area,fni
471
472
473
474 INTEGER I, IL, LOLD, JJ, NN, J3,J2, J1, IG, I3, I2, I1, K, IFQ, MFROT, NISKYL
476 . h(4), xx1(4), xx2(4), xx3(4),
477 . fxi(mvsiz), fyi(mvsiz),
478 . fzi(mvsiz), fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz), fy1(mvsiz), fy2(mvsiz),
479 . fy3(mvsiz), fy4(mvsiz), fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
480 . ss0, tt0, xc, econvt,
alpha, alphi,
481 . yc, zc, xc0, yc0, zc0, sp, sm, tp, tm, ansx, ansy, ansz, fmax,
482 . stf, fti, fn, tn1, tn2, tn3, tn, dtm, xmu, vx,vy,vz,vv,v2,p,
483 . vv1,vv2,dmu,aa
484 INTEGER IRTLO_I(MVSIZ)
486
487
488
489
490 IF (inconv/=1) THEN
491 DO i=lft,llt
492 il=i+nft
493 irtlo_i(i)=irtlo(il)
494 fric0_i(1,i)=fric0(1,il)
495 fric0_i(2,i)=fric0(2,il)
496 fric0_i(3,i)=fric0(3,il)
497 ENDDO
498 ENDIF
499 DO 300 i=lft,llt
500 fxi(i)=zero
501 fyi(i)=zero
502 fzi(i)=zero
503 il=i+nft
504 econvt = zero
505 IF(xface(i)==zero) THEN
506
507
508
509 irtlo(il)=0
510 fric0(1,il)=zero
511 fric0(2,il)=zero
512 fric0(3,il)=zero
513 ELSE
514
515 lold=iabs(irtlo(il))
516 IF(lold==0)THEN
517
518
519
520 irtlo(il)=irtl(il)*xface(i)
521 cst(1,il)=ssc(i)
522 cst(2,il)=ttc(i)
523 ELSE
524
525
526
527 ss0=cst(1,il)
528 tt0=cst(2,il)
529 fxi(i)=fric0(1,il)
530 fyi(i)=fric0(2,il)
531 fzi(i)=fric0(3,il)
532
533 xc=xp(i)
534 yc=yp(i)
535 zc=zp(i)
536 DO jj=1,4
537 nn=msr(irect(jj,lold))
538 xx1(jj)=x(1,nn)
539 xx2(jj)=x(2,nn)
540 xx3(jj)=x(3,nn)
541 ENDDO
542 xc0=zero
543 yc0=zero
544 zc0=zero
545 sp=one+ss0
546 sm=one-ss0
547 tp=fourth*(one + tt0)
548 tm=fourth*(one - tt0)
549 h(1)=tm*sm
550 h(2)=tm*sp
551 h(3)=tp*sp
552 h(4)=tp*sm
553 DO jj=1,4
554 xc0=xc0+h(jj)*xx1(jj)
555 yc0=yc0+h(jj)*xx2(jj)
556 zc0=zc0+h(jj)*xx3(jj)
557 ENDDO
558 ansx= (xc-xc0)
559 ansy= (yc-yc0)
560 ansz= (zc-zc0)
561 econvt = econvt + half*(fxi(i)*ansx+fyi(i)*ansy+fzi(i)*ansz)
562 stf=em01*stif(i)
563 fxi(i)=fxi(i) + ansx*stf
564 fyi(i)=fyi(i) + ansy*stf
565 fzi(i)=fzi(i) + ansz*stf
566
567 IF (codvers>=44) THEN
568 ifq = ipari(31)
569 IF (ifq>0) THEN
570 IF (ifq==3) freq =
max(one,freq*dt12)
573 k = 3*(il-1)
574 IF (fni(i)==0) THEN
575 ftsav(k+1) = zero
576 ftsav(k+2) = zero
577 ftsav(k+3) = zero
578 ELSE
579 fxi(i)=
alpha*fxi(i) + alphi*ftsav(k+1)
580 fyi(i)=
alpha*fyi(i) + alphi*ftsav(k+2)
581 fzi(i)=
alpha*fzi(i) + alphi*ftsav(k+3)
582 IF (inconv==1) THEN
583 ftsav(k+1) = fxi(i)
584 ftsav(k+2) = fyi(i)
585 ftsav(k+3) = fzi(i)
586 ENDIF
587 ENDIF
588 ENDIF
589 ENDIF
590 fti=sqrt(fxi(i)*fxi(i)+fyi(i)*fyi(i)+fzi(i)*fzi(i))
591
592 fn =fxi(i)*n1(i)+fyi(i)*n2(i)+fzi(i)*n3(i)
593 tn1=fxi(i)-n1(i)*fn
594 tn2=fyi(i)-n2(i)*fn
595 tn3=fzi(i)-n3(i)*fn
596 tn =sqrt(tn1*tn1+tn2*tn2+tn3*tn3)
597 IF(tn/=zero)THEN
598 tn1=tn1/tn
599 tn2=tn2/tn
600 tn3=tn3/tn
601 ELSE
602 tn3=zero
603 tn=sqrt(n1(i)*n1(i)+n2(i)*n2(i))
604 IF(tn/=zero)THEN
605 tn2=-n1(i)/tn
606 tn1= n2(i)/tn
607 ELSE
608 tn2=zero
609 tn1=one
610 ENDIF
611 ENDIF
612
614 xmu = fric
615 IF(codvers>=44) THEN
616
617
618
619 mfrot = ipari(30)
620 IF(mfrot>=1) THEN
621 ig=nsv(il)
622 vx = v(1,ig) - h(1)*v(1,ix1(i)) - h(2)*v(1,ix2(i))
623 . - h(3)*v(1,ix3(i)) - h(4)*v(1,ix4(i))
624 vy = v(2,ig) - h(1)*v(2,ix1(i)) - h(2)*v(2,ix2(i))
625 . - h(3)*v(2,ix3(i)) - h(4)*v(2,ix4(i))
626 vz = v(3,ig) - h(1)*v(3,ix1(i)) - h(2)*v(3,ix2(i))
627 . - h(3)*v(3,ix3(i)) - h(4)*v(3,ix4(i))
628
629 aa = n1(i)*vx + n2(i)*vy + n3(i)*vz
630 vx = vx - n1(i)*aa
631 vy = vy - n2(i)*aa
632 vz = vz - n3(i)*aa
633 v2 = vx*vx + vy*vy + vz*vz
634 vv = sqrt(
max(em30, v2))
635 ENDIF
636 IF(mfrot==1)THEN
637 xmu = fric + (frot_p(1) + frot_p(4)*p ) * p
638 . +(frot_p(2) + frot_p(3)*p) * vv + frot_p(5)*v2
639 ELSEIF(mfrot==2)THEN
640
641 xmu = fric
642 . + frot_p(1)*exp(frot_p(2)*vv)*p**2
643 . + frot_p(3)*exp(frot_p(4)*vv)*p
644 . + frot_p(5)*exp(frot_p(6)*vv)
645 ELSEIF(mfrot==3)THEN
646
647 IF(vv>=0.AND.vv<=frot_p(5)) THEN
648 dmu = frot_p(3)-frot_p(1)
649 vv1 = vv / frot_p(5)
650 xmu = frot_p(1)+ dmu*vv1*(two-vv1)
651 ELSEIF(vv>frot_p(5).AND.vv<frot_p(6)) THEN
652 dmu = frot_p(4)-frot_p(3)
653 vv1 = (vv - frot_p(5))/(frot_p(6)-frot_p(5))
654 xmu = frot_p(3)+ dmu * (three-two*vv1)*vv1**2
655 ELSE
656 dmu = frot_p(2)-frot_p(4)
657 vv2 = (vv - frot_p(6))**2
658 xmu = frot_p(2) - dmu / (one + dmu*vv2)
659 ENDIF
660 ENDIF
661 fmax = -
min(xmu*fni(i),zero)
662 IF(fti>fmax)THEN
663
664 fxi(i)=tn1*fmax
665 fyi(i)=tn2*fmax
666 fzi(i)=tn3*fmax
667 fric0(1,il)=fxi(i)
668 fric0(2,il)=fyi(i)
669 fric0(3,il)=fzi(i)
670 irtlo(il)=irtl(il)*xface(i)
671 cst(1,il)=ssc(i)
672 cst(2,il)=ttc(i)
673 ELSE
674
675 fxi(i)=tn1*fti
676 fyi(i)=tn2*fti
677 fzi(i)=tn3*fti
678 ENDIF
679 ELSE
680
681
682
683 ig=nsv(il)
684 vx = v(1,ig) - h(1)*v(1,ix1(i)) - h(2)*v(1,ix2(i))
685 . - h(3)*v(1,ix3(i)) - h(4)*v(1,ix4(i))
686 vy = v(2,ig) - h(1)*v(2,ix1(i)) - h(2)*v(2,ix2(i))
687 . - h(3)*v(2,ix3(i)) - h(4)*v(2,ix4(i))
688 vz = v(3,ig) - h(1)*v(3,ix1(i)) - h(2)*v(3,ix2(i))
689 . - h(3)*v(3,ix3(i)) - h(4)*v(3,ix4(i))
690 vv = vx*vx + vy*vy + vz*vz
691 xmu = fric + ( cf(1) + cf(4)*p ) * p
692 . + ( cf(2) + cf(3)*p ) * sqrt(vv) + cf(5)*vv
693 fmax= -
min(xmu*fni(i),zero)
694 IF(fti>fmax)THEN
695
696
697
698 fxi(i)=tn1*fmax
699 fyi(i)=tn2*fmax
700 fzi(i)=tn3*fmax
701 fric0(1,il)=fxi(i)
702 fric0(2,il)=fyi(i)
703 fric0(3,il)=fzi(i)
704 irtlo(il)=irtl(il)*xface(i)
705 cst(1,il)=ssc(i)
706 cst(2,il)=ttc(i)
707 ELSE
708
709
710
711 fxi(i)=tn1*fti
712 fyi(i)=tn2*fti
713 fzi(i)=tn3*fti
714 ENDIF
715 ENDIF
716 econvt = econvt + half*(fxi(i)*ansx+fyi(i)*ansy+fzi(i)*ansz)
717 ENDIF
718 ENDIF
719
720 300 CONTINUE
721
722 IF (inconv==1) THEN
723
724 econtv = econtv + econvt
725
726 fsav(27) = fsav(27) + econvt
727
728
729
730 dtm=imast*dt12
731 DO i=lft,llt
732 fsav(4)=fsav(4)+fxi(i)*dtm
733 fsav(5)=fsav(5)+fyi(i)*dtm
734 fsav(6)=fsav(6)+fzi(i)*dtm
735 ENDDO
736 END IF
737
738 DO i=lft,llt
739 fx1(i)=fxi(i)*h1(i)
740 fy1(i)=fyi(i)*h1(i)
741 fz1(i)=fzi(i)*h1(i)
742
743 fx2(i)=fxi(i)*h2(i)
744 fy2(i)=fyi(i)*h2(i)
745 fz2(i)=fzi(i)*h2(i)
746
747 fx3(i)=fxi(i)*h3(i)
748 fy3(i)=fyi(i)*h3(i)
749 fz3(i)=fzi(i)*h3(i)
750
751 fx4(i)=fxi(i)*h4(i)
752 fy4(i)=fyi(i)*h4(i)
753 fz4(i)=fzi(i)*h4(i)
754 ENDDO
755
756 IF(iparit==0)THEN
757
758
759
760 DO i=lft,llt
761 il=i+nft
762 j3=3*ix1(i)
763 j2=j3-1
764 j1=j2-1
765 e(j1)=e(j1)+fx1(i)
766 e(j2)=e(j2)+fy1(i)
767 e(j3)=e(j3)+fz1(i)
768
769 j3=3*ix2(i)
770 j2=j3-1
771 j1=j2-1
772 e(j1)=e(j1)+fx2(i)
773 e(j2)=e(j2)+fy2(i)
774 e(j3)=e(j3)+fz2(i)
775
776 j3=3*ix3(i)
777 j2=j3-1
778 j1=j2-1
779 e(j1)=e(j1)+fx3(i)
780 e(j2)=e(j2)+fy3(i)
781 e(j3)=e(j3)+fz3(i)
782
783 j3=3*ix4(i)
784 j2=j3-1
785 j1=j2-1
786 e(j1)=e(j1)+fx4(i)
787 e(j2)=e(j2)+fy4(i)
788 e(j3)=e(j3)+fz4(i)
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)=
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>=output%TANIM .AND. tt<=output%TANIM_STOP).OR.tt>=toutp.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>=output%TANIM .AND. tt<=output%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)