35
36
37
38#include "implicit_f.inc"
39
40
41
42#include "mvsiz_p.inc"
43
44
45
46#include "com08_c.inc"
47
48
49
50 INTEGER, INTENT(IN) :: NEL
51 INTEGER NC1(MVSIZ),NC2(MVSIZ)
53 . v(3,*),offg(mvsiz),off(mvsiz),eps(mvsiz),al(mvsiz),
54 . vx1(mvsiz),vx2(mvsiz),vy1(mvsiz),vy2(mvsiz),vz1(mvsiz),
55 . vz2(mvsiz),ex(mvsiz),ey(mvsiz),ez(mvsiz),x1(mvsiz),x2(mvsiz),
56 . y1(mvsiz),y2(mvsiz),z1(mvsiz),z2(mvsiz)
57
58
59
60 INTEGER I,J
62 . v1(mvsiz),v2(mvsiz),ratio(mvsiz),off_l
63
64 DO i=1,nel
65 vx1(i)=v(1,nc1(i))
66 vy1(i)=v(2,nc1(i))
67 vz1(i)=v(3,nc1(i))
68 vx2(i)=v(1,nc2(i))
69 vy2(i)=v(2,nc2(i))
70 vz2(i)=v(3,nc2(i))
71 ENDDO
72
73 DO i=1,nel
74 ex(i)=x2(i)-x1(i)
75 ey(i)=y2(i)-y1(i)
76 ez(i)=z2(i)-z1(i)
77 ENDDO
78
79 DO i=1,nel
80 al(i)=sqrt(ex(i)*ex(i)+ey(i)*ey(i)+ez(i)*ez(i))
81 ENDDO
82
83 DO i=1,nel
84 IF (al(i) /= zero) THEN
85 ex(i) = ex(i)/al(i)
86 ey(i) = ey(i)/al(i)
87 ez(i) = ez(i)/al(i)
88 ELSE
89 ex(i) = zero
90 ey(i) = zero
91 ez(i) = zero
92 ENDIF
93 ENDDO
94
95 DO i=1,nel
96 v1(i)=ex(i)*vx1(i)+ey(i)*vy1(i)+ez(i)*vz1(i)
97 v2(i)=ex(i)*vx2(i)+ey(i)*vy2(i)+ez(i)*vz2(i)
98 ENDDO
99
100 DO i=1,nel
101 IF (al(i) /= zero) THEN
102 eps(i)= (v2(i)-v1(i))/al(i)
103 ratio(i) = ( (vx2(i)-vx1(i))*(vx2(i)-vx1(i))
104 . +(vy2(i)-vy1(i))*(vy2(i)-vy1(i))
105 . +(vz2(i)-vz1(i))*(vz2(i)-vz1(i)) )
106 . /(al(i)*al(i))
107 ELSE
108 eps(i) = zero
109 ratio(i) = zero
110 ENDIF
111 ENDDO
112
113 DO i=1,nel
114 eps(i)= eps(i) + half*dt1*( eps(i)*eps(i) - ratio(i) )
115 ENDDO
116
117 off_l = zero
118 DO i=1,nel
119 off(i) =
min(one,abs(offg(i)))
120 off(i) =
max(zero,off(i))
121 off_l =
min(off_l,offg(i))
122 ENDDO
123
124 IF (off_l < zero) THEN
125 DO i=1,nel
126 IF (offg(i) < zero) THEN
127 eps(i)=zero
128 ENDIF
129 ENDDO
130 ENDIF
131
132 RETURN