45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57
58
59
60 INTEGER, INTENT(IN) ::
61 INTEGER, INTENT(IN) :: JHBE
63 . off(*),det(*),
64 . x1(*), x2(*), x3(*), x4(*), x5(*), x6(*), x7(*), x8(*),
65 . y1(*), y2(*), y3(*), y4(*), y5(*), y6(*), y7(*), y8(*),
66 . z1(*), z2(*), z3(*), z4(*), z5(*), z6(*), z7(*), z8(*),
67 . px1(*), px2(*), px3(*), px4(*),
68 . py1(*), py2(*), py3(*), py4(*),
69 . pz1(*), pz2(*), pz3(*), pz4(*),
70 . px1h1(*), px1h2(*), px1h3(*),
71 . px2h1(*), px2h2(*), px2h3(*),
72 . px3h1(*), px3h2(*), px3h3(*),
73 . px4h1(*), px4h2(*), px4h3(*)
74
75
76
77 INTEGER NGL(*), I
79 . dett(mvsiz) , aj7(mvsiz), aj8(mvsiz) , aj9(mvsiz),
80 . aji1(mvsiz), aji2(mvsiz), aji3(mvsiz),
81 . aji4(mvsiz), aji5(mvsiz), aji6(mvsiz),
82 . aji7(mvsiz), aji8(mvsiz), aji9(mvsiz),
83 . x17(mvsiz) , x28(mvsiz) , x35(mvsiz) , x46(mvsiz),
84 . y17(mvsiz) , y28(mvsiz) , y35(mvsiz) , y46(mvsiz),
85 . z17(mvsiz) , z28(mvsiz) , z35(mvsiz) , z46(mvsiz),
86 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
87 . aj12(mvsiz), aj45(mvsiz), aj78(mvsiz),
88 . a17(mvsiz) , a28(mvsiz) ,
89 . b17(mvsiz) , b28(mvsiz) ,
90 . c17(mvsiz) , c28(mvsiz) , aj4(mvsiz),
91 . aj5(mvsiz) , aj6(mvsiz) , aj1(mvsiz),
92 . aj2(mvsiz) , aj3(mvsiz) ,
93 . hx,hy,hz
94
95
96 DO i=1,nel
97 x17(i)=x7(i)-x1(i)
98 x28(i)=x8(i)-x2(i)
99 x35(i)=x5(i)-x3(i)
100 x46(i)=x6(i)-x4(i)
101 y17(i)=y7(i)-y1(i)
102 y28(i)=y8(i)-y2(i)
103 y35(i)=y5(i)-y3(i)
104 y46(i)=y6(i)-y4(i)
105 z17(i)=z7(i)-z1(i)
106 z28(i)=z8(i)-z2(i)
107 z35(i)=z5(i)-z3(i)
108 z46(i)=z6(i)-z4(i)
109 ENDDO
110
111 DO i=1,nel
112 aj1(i)=x17(i)+x28(i)-x35(i)-x46(i)
113 aj2(i)=y17(i)+y28(i)-y35(i)-y46(i)
114 aj3(i)=z17(i)+z28(i)-z35(i)-z46(i)
115 a17(i)=x17(i)+x46(i)
116 a28(i)=x28(i)+x35(i)
117 b17(i)=y17(i)+y46(i)
118 b28(i)=y28(i)+y35(i)
119 c17(i)=z17(i)+z46(i)
120 c28(i)=z28(i)+z35(i)
121 ENDDO
122
123 DO i=1,nel
124 aj4(i)=a17(i)+a28(i)
125 aj5(i)=b17(i)+b28(i)
126 aj6(i)=c17(i)+c28(i)
127 aj7(i)=a17(i)-a28(i)
128 aj8(i)=b17(i)-b28(i)
129 aj9(i)=c17(i)-c28(i)
130 ENDDO
131
132
133
134 DO i=1,nel
135 jac_59_68(i)=aj5(i)*aj9(i)-aj6(i)*aj8(i)
136 jac_67_49(i)=aj6(i)*aj7(i)-aj4(i)*aj9(i)
137 jac_48_57(i)=aj4(i)*aj8(i)-aj5(i)*aj7(i)
138 ENDDO
139
140 DO i=1,nel
141 det(i)=one_over_64*(aj1(i)*jac_59_68(i)+aj2(i)*jac_67_49(i)+aj3(i)*jac_48_57(i))
142 ENDDO
143
145 1 off, det, ngl, nel)
146
147 DO i=1,nel
148 dett(i)=one_over_64/det(i)
149 ENDDO
150
151
152
153
154
155 DO i=1,nel
156 aji1(i)=dett(i)*jac_59_68(i)
157 aji4(i)=dett(i)*jac_67_49(i)
158 aji7(i)=dett(i)*jac_48_57(i)
159 aji2(i)=dett(i)*(-aj2(i)*aj9(i)+aj3(i)*aj8(i))
160 aji5(i)=dett(i)*( aj1(i)*aj9(i)-aj3(i)*aj7(i))
161 aji8(i)=dett(i)*(-aj1(i)*aj8(i)+aj2(i)*aj7(i))
162 aji3(i)=dett(i)*( aj2(i)*aj6(i)-aj3(i)*aj5(i))
163 aji6(i)=dett(i)*(-aj1(i)*aj6(i)+aj3(i)*aj4(i))
164 aji9(i)=dett(i)*( aj1(i)*aj5(i)-aj2(i)*aj4(i))
165 ENDDO
166
167 DO i=1,nel
168 aj12(i)=aji1(i)-aji2(i)
169 aj45(i)=aji4(i)-aji5(i)
170 aj78(i)=aji7(i)-aji8(i)
171
172 px3(i)= aj12(i)+aji3(i)
173 py3(i)= aj45(i)+aji6(i)
174 pz3(i)= aj78(i)+aji9(i)
175 px4(i)= aj12(i)-aji3(i)
176 py4(i)= aj45(i)-aji6(i)
177 pz4(i)= aj78(i)-aji9(i)
178
179 aj12(i)=aji1(i)+aji2(i)
180 aj45(i)=aji4(i)+aji5(i)
181 aj78(i)=aji7(i)+aji8(i)
182
183 px1(i)=-aj12(i)-aji3(i)
184 py1(i)=-aj45(i)-aji6(i)
185 pz1(i)=-aj78(i)-aji9(i)
186 px2(i)=-aj12(i)+aji3(i)
187 py2(i)=-aj45(i)+aji6(i)
188 pz2(i)=-aj78(i)+aji9(i)
189 ENDDO
190
191 IF(jhbe/=0)THEN
192 DO i=1,nel
193
194 hx=(x1(i)-x2(i)+x3(i)-x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
195 hy=(y1(i)-y2(i)+y3(i)-y4(i)+y5(i)-y6(i)+y7(i)-y8(i))
196 hz=(z1(i)-z2(i)+z3(i)-z4(i)+z5(i)-z6(i)+z7(i)-z8(i))
197 px1h1(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
198 px2h1(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
199 px3h1(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
200 px4h1(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
201 ENDDO
202
203 DO i=1,nel
204 hx=(x1(i)+x2(i)-x3(i)-x4(i)-x5(i)-x6(i)+x7(i)+x8(i))
205 hy=(y1(i)+y2(i)-y3(i)-y4(i)-y5(i)-y6(i)+y7(i)+y8(i))
206 hz=(z1(i)+z2(i)-z3(i)-z4(i)-z5(i)-z6(i)+z7(i)+z8(i))
207 px1h2(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
208 px2h2(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
209 px3h2(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
210 px4h2(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
211 ENDDO
212
213
214 DO i=1,nel
215 hx=(x1(i)-x2(i)-x3(i)+x4(i)-x5(i)+x6(i)+x7(i)-x8(i))
216 hy=(y1(i)-y2(i)-y3(i)+y4(i)-y5(i)+y6(i)+y7(i)-y8(i))
217 hz=(z1(i)-z2(i)-z3(i)+z4(i)-z5(i)+z6(i)+z7(i)-z8(i))
218 px1h3(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
219 px2h3(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
220 px3h3(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
221 px4h3(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
222 ENDDO
223
224 ENDIF
225 RETURN
226
subroutine schkjab3(off, det, ngl, nel)