37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "mvsiz_p.inc"
45
46
47
48#include "com08_c.inc"
49
50
51
52 INTEGER, INTENT(IN) :: NEL
53 INTEGER NC1(MVSIZ), NC2(MVSIZ), NC3(MVSIZ)
55 . rloc(3,*),r(3,*),al(*),x1(mvsiz),x2(mvsiz),x3(mvsiz),
56 . y1(mvsiz),y2(mvsiz),y3(mvsiz),z1(mvsiz),z2(mvsiz),
57 . z3(mvsiz),rx1g(mvsiz),rx2g(mvsiz),ry1g(mvsiz),ry2g(mvsiz),
58 . rz1g(mvsiz),rz2g(mvsiz),rx1(mvsiz),rx2(mvsiz),ry1(mvsiz),ry2(mvsiz),
59 . rz1(mvsiz),rz2(mvsiz),e1x(mvsiz),e1y(mvsiz),e1z(mvsiz),e2x(mvsiz),
60 . e2y(mvsiz),e2z(mvsiz),e3x(mvsiz),e3y(mvsiz),e3z(mvsiz)
61
62
63
64 INTEGER I, J
66 . sum(mvsiz),
67 . theta, sum2(mvsiz), sum3(mvsiz), cost(mvsiz),
68 . sint(mvsiz), r1phi, r2phi, thetaphi, sum2phi, sum3phi, cophi,
69 . siphi, sinphi, cosphi, sumphi
70
71 DO i=1,nel
72 rx1g(i)=r(1,nc1(i))
73 ry1g(i)=r(2,nc1(i))
74 rz1g(i)=r(3,nc1(i))
75 rx2g(i)=r(1,nc2(i))
76 ry2g(i)=r(2,nc2(i))
77 rz2g(i)=r(3,nc2(i))
78 ENDDO
79
80 DO i=1,nel
81 e2x(i)=rloc(1,i)
82 e2y(i)=rloc(2,i)
83 e2z(i)=rloc(3,i)
84 ENDDO
85
86 DO i=1,nel
87 e1x(i)=x2(i)-x1(i)
88 e1y(i)=y2(i)-y1(i)
89 e1z(i)=z2(i)-z1(i)
90 ENDDO
91
92 DO i=1,nel
93 al(i)=sqrt(e1x(i)**2+e1y(i)**2+e1z(i)**2)
94 ENDDO
95
96 DO i=1,nel
97 e1x(i)=e1x(i)/al(i)
98 e1y(i)=e1y(i)/al(i)
99 e1z(i)=e1z(i)/al(i)
100 ENDDO
101
102 DO i=1,nel
103 e3x(i)=e1y(i)*e2z(i)-e1z(i)*e2y(i)
104 e3y(i)=e1z(i)*e2x(i)-e1x(i)*e2z(i)
105 e3z(i)=e1x(i)*e2y(i)-e1y(i)*e2x(i)
106 ENDDO
107
108 DO i=1,nel
109 e2x(i)=e3y(i)*e1z(i)-e3z(i)*e1y(i)
110 e2y(i)=e3z(i)*e1x(i)-e3x(i)*e1z(i)
111 e2z(i)=e3x(i)*e1y(i)-e3y(i)*e1x(i)
112 ENDDO
113
114
115
116 DO i=1,nel
117 rx1(i)=e1x(i)*rx1g(i)+e1y(i)*ry1g(i)+e1z(i)*rz1g(i)
118 rx2(i)=e1x(i)*rx2g(i)+e1y(i)*ry2g(i)+e1z(i)*rz2g(i)
119 theta=(rx1(i)+rx2(i))/two*dt1
120 sum2(i)=sqrt(e2x(i)**2+e2y(i)**2+e2z(i)**2)
121 sum3(i)=sqrt(e3x(i)**2+e3y(i)**2+e3z(i)**2)
122 cost(i)=cos(theta)/sum2(i)
123 sint(i) =sin(theta)/sum3(i)
124 ENDDO
125
126
127
128 DO i=1,nel
129 e2x(i)=e2x(i)*cost(i)+e3x(i)*sint(i)
130 e2y(i)=e2y(i)*cost(i)+e3y(i)*sint(i)
131 e2z(i)=e2z(i)*cost(i)+e3z(i)*sint(i)
132 ENDDO
133
134 DO i=1,nel
135 sum(i)=sqrt(e2x(i)**2+e2y(i)**2+e2z(i)**2)
136 ENDDO
137
138 DO i=1,nel
139 e2x(i)=e2x(i)/sum(i)
140 e2y(i)=e2y(i)/sum(i)
141 e2z(i)=e2z(i)/sum(i)
142 ENDDO
143
144 DO i=1,nel
145 e3x(i)=e1y(i)*e2z(i)-e1z(i)*e2y(i)
146 e3y(i)=e1z(i)*e2x(i)-e1x(i)*e2z(i)
147 e3z(i)=e1x(i)*e2y(i)-e1y(i)*e2x(i)
148 ENDDO
149
150 DO i=1,nel
151 sum(i)=sqrt(e3x(i)**2+e3y(i)**2+e3z(i)**2)
152 e3x(i)=e3x(i)/sum(i)
153 e3y(i)=e3y(i)/sum(i)
154 e3z(i)=e3z(i)/sum(i)
155 ENDDO
156
157 DO i=1,nel
158 rloc(1,i)=e2x(i)
159 rloc(2,i)=e2y(i)
160 rloc(3,i)=e2z(i)
161 ENDDO
162
163 RETURN