35
36
37
38#include "implicit_f.inc"
39
40
41
42#include "vect01_c.inc"
43#include "param_c.inc"
44
45
46
47 INTEGER NEL,NUVAR,IREP
48
49 my_real uparam(*),uvar(nel,*),aldt(nel),tan_phi(nel),dir1(nel,2),dir2(nel,2)
50 my_real,
DIMENSION(NEL),
INTENT(IN) :: e1x,e2x,e3x,e1y,e2y,e3y,e1z,e2z,e3z,
51 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
52
53
54
55 INTEGER I
57 . r,s,d1,d2,d11,d12,d21,d22,u1x,u1y,u2x,u2y,det,w1x,w2x,w1y,w2y,
58 . phi,csp,snp,tana,r1,r2,s1,s2,g0,aa,bb,suma,v1,v2,v3,wr,ws
59 my_real ,
DIMENSION(NEL) :: e11,e12,e13,e21,e22,e23
60
61
62 IF (ity == 3) THEN
63
64 DO i=1,nel
65 e11(i)= x2(i)+x3(i)-x1(i)-x4(i)
66 e12(i)= y2(i)+y3(i)-y1(i)-y4(i)
67 e13(i)= z2(i)+z3(i)-z1(i)-z4(i)
68 e21(i)= x3(i)+x4(i)-x1(i)-x2(i)
69 e22(i)= y3(i)+y4(i)-y1(i)-y2(i)
70 e23(i)= z3(i)+z4(i)-z1(i)-z2(i)
71 ENDDO
72 ELSEIF (ity == 7) THEN
73
74 DO i=1,nel
75 e11(i)= x2(i)-x1(i)
76 e12(i)= y2(i)-y1(i)
77 e13(i)= z2(i)-z1(i)
78 e21(i)= x3(i)-x1(i)
79 e22(i)= y3(i)-y1(i)
80 e23(i)= z3(i)-z1(i)
81 ENDDO
82 ENDIF
83
84 g0 = uparam(13)
85 DO i=1,nel
86 aa = dir1(i,1)
87 bb = dir1(i,2)
88 v1 = aa*e11(i) + bb*e21(i)
89 v2 = aa*e12(i) + bb*e22(i)
90 v3 = aa*e13(i) + bb*e23(i)
91 wr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
92 ws = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
93 suma =
max( sqrt(wr*wr + ws*ws), em20)
94 r1 = wr/suma
95 s1 = ws/suma
96 aa = dir2(i,1)
97 bb = dir2(i,2)
98 v1 = aa*e11(i) + bb*e21(i)
99 v2 = aa*e12(i) + bb*e22(i)
100 v3 = aa
101 wr = v1*e1x(i)+ v2*e1y(i) + v3*e1z(i)
102 ws = v1*e2x(i)+ v2*e2y(i) + v3*e2z(i)
103 suma =
max( sqrt(wr*wr + ws*ws), em20)
104 r2 = wr/suma
105 s2 = ws/suma
106 tana = (r1*r2 + s1*s2) / (r1*s2 - r2*s1)
107
108 uvar(i,1:nuvar) = zero
109 uvar(i,6) = tana
110 tan_phi(i) = tana
111 uvar(i,10) = tana*g0
112 uvar(i,14) = aldt(i)
113 uvar(i,40) = one
114 ENDDO
115
116 RETURN