48
49
50
51#include "implicit_f.inc"
52#include "comlock.inc"
53
54
55
56#include "mvsiz_p.inc"
57
58
59
60 INTEGER, INTENT(IN) :: JLAG
61 INTEGER,INTENT(IN) :: NGL(*),NEL,JHBE,ISMSTR
62 double precision
63 . x1(*), x2(*), x3(*), x4(*), x5(*), x6(*), x7(*), x8(*),
64 . y1(*), y2(*), y3(*), y4(*), y5(*), y6(*), y7(*), y8(*),
65 . z1(*), z2(*), z3(*), z4(*), z5(*), z6(*), z7(*), z8(*),
66 . sav(nel,21) ,voldp(*)
68 . off(*), det(*),
69 . px1(*), px2(*), px3(*), px4(*),
70 . py1(*), py2(*), py3(*), py4(*),
71 . pz1(*), pz2(*), pz3(*), pz4(*),
72 . px1h1(*), px1h2(*), px1h3(*),
73 . px2h1(*), px2h2(*), px2h3(*),
74 . px3h1(*), px3h2(*), px3h3(*),
75 . px4h1(*), px4h2(*), px4h3(*),
76 . jac1(*), jac2(*), jac3(*),
77 . jac4(*), jac5(*), jac6(*), offg(*)
78
79
80
81 INTEGER I, J, NNEGA, INDEX(MVSIZ)
82
84 . dett , jac7(mvsiz), jac8(mvsiz) , jac9(mvsiz),
85 . jaci1, jaci2, jaci3,
86 . jaci4, jaci5, jaci6,
87 . jaci7, jaci8, jaci9,
88 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
89 . x_17_46, x_28_35,
90 . y_17_46, y_28_35,
91 . z_17_46, z_28_35,
92 . jac12, jac45, jac78,
93 . hx, hy, hz
94 double precision
95 . x17 , x28 , x35 , x46,
96 . y17 , y28 , y35 , y46,
97 . z17 , z28 , z35 , z46,
98 . a17_46 , a28_35,
99 . b17_46 , b28_35 ,
100 . c17_46 , c28_35
101
102
103
104 DO i=1,nel
105 x17=x7(i)-x1(i)
106 x28=x8(i)-x2(i)
107 x35=x5(i)-x3(i)
108 x46=x6(i)-x4(i)
109
110 y17=y7(i)-y1(i)
111 y28=y8(i)-y2(i)
112 y35=y5(i)-y3(i)
113 y46=y6(i)-y4(i)
114
115 z17=z7(i)-z1(i)
116 z28=z8(i)-z2(i)
117 z35=z5(i)-z3(i)
118 z46=z6(i)-z4(i)
119
120 jac1(i)=x17+x28-x35-x46
121 jac2(i)=y17+y28-y35-y46
122 jac3(i)=z17+z28-z35-z46
123
124 x_17_46=x17+x46
125 x_28_35=x28+x35
126 y_17_46=y17+y46
127 y_28_35=y28+y35
128 z_17_46=z17+z46
129 z_28_35=z28+z35
130
131 jac4(i)=x_17_46+x_28_35
132 jac5(i)=y_17_46+y_28_35
133 jac6(i)=z_17_46+z_28_35
134
135 jac7(i)=x_17_46-x_28_35
136 jac8(i)=y_17_46-y_28_35
137 jac9(i)=z_17_46-z_28_35
138 ENDDO
139
140 DO i=1,nel
141 jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
142 jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
143 jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
144 ENDDO
145
146 DO i=1,nel
147 voldp(i)=one_over_64*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
148 det(i)=voldp(i)
149 ENDDO
150
151
153 1 off, det, ngl, offg,
154 2 nnega, index, nel, ismstr,
155 3 jlag)
156 IF (nnega>0) THEN
157 IF (ismstr==10.OR.ismstr==12) THEN
158#include "vectorize.inc"
159 DO j=1,nnega
160 i = index(j)
161 x1(i)=sav(i,1)
162 y1(i)=sav(i,8)
163 z1(i)=sav(i,15)
164 x2(i)=sav(i,2)
165 y2(i)=sav(i,9)
166 z2(i)=sav(i,16)
167 x3(i)=sav(i,3)
168 y3(i)=sav(i,10)
169 z3(i)=sav(i,17)
170 x4(i)=sav(i,4)
171 y4(i)=sav(i,11)
172 z4(i)=sav(i,18)
173 x5(i)=sav(i,5)
174 y5(i)=sav(i,12)
175 z5(i)=sav(i,19)
176 x6(i)=sav(i,6)
177 y6(i)=sav(i,13)
178 z6(i)=sav(i,20)
179 x7(i)=sav(i,7)
180 y7(i)=sav(i,14)
181 z7(i)=sav(i,21)
182 x8(i)=zero
183 y8(i)=zero
184 z8(i)=zero
185 ENDDO
186 ELSE
187#include "vectorize.inc"
188 DO j=1,nnega
189 i = index(j)
190 x1(i)=sav(i,1)
191 y1(i)=sav(i,2)
192 z1(i)=sav(i,3)
193 x2(i)=sav(i,4)
194 y2(i)=sav(i,5)
195 z2(i)=sav(i,6)
196 x3(i)=sav(i,7)
197 y3(i)=sav(i,8)
198 z3(i)=sav(i,9)
199 x4(i)=sav(i,10)
200 y4(i)=sav(i,11)
201 z4(i)=sav(i,12)
202 x5(i)=sav(i,13)
203 y5(i)=sav(i,14)
204 z5(i)=sav(i,15)
205 x6(i)=sav(i,16)
206 y6(i)=sav(i,17)
207 z6(i)=sav(i,18)
208 x7(i)=sav(i,19)
209 y7(i)=sav(i,20)
210 z7(i)=sav(i,21)
211 x8(i)=zero
212 y8(i)=zero
213 z8(i)=zero
214 ENDDO
215 END IF
216#include "vectorize.inc"
217 DO j=1,nnega
218 i = index(j)
219
220 x17=x7(i)-x1(i)
221 x28=x8(i)-x2(i)
222 x35=x5(i)-x3(i)
223 x46=x6(i)-x4(i)
224
225 y17=y7(i)-y1(i)
226 y28=y8(i)-y2(i)
227 y35=y5(i)-y3(i)
228 y46=y6(i)-y4(i)
229
230 z17=z7(i)-z1(i)
231 z28=z8(i)-z2(i)
232 z35=z5(i)-z3(i)
233 z46=z6(i)-z4(i)
234
235 jac1(i)=x17+x28-x35-x46
236 jac2(i)=y17+y28-y35-y46
237 jac3(i)=z17+z28-z35-z46
238
239 a17_46=x17+x46
240 a28_35=x28+x35
241 b17_46=y17+y46
242 b28_35=y28+y35
243 c17_46=z17+z46
244 c28_35=z28+z35
245
246 jac4(i)=a17_46+a28_35
247 jac5(i)=b17_46+b28_35
248 jac6(i)=c17_46+c28_35
249 jac7(i)=a17_46-a28_35
250 jac8(i)=b17_46-b28_35
251 jac9(i)=c17_46-c28_35
252
253 jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
254 jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
255 jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
256
257 voldp(i)=one_over_64*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
258 det(i)=voldp(i)
259
260 offg(i) = two
261 ENDDO
262 END IF
263
264
265 DO i=1,nel
266 dett=one_over_64/det(i)
267 jaci1=dett*jac_59_68(i)
268 jaci4=dett*jac_67_49(i)
269 jaci7=dett*jac_48_57(i)
270 jaci2=dett*(jac3(i)*jac8(i)-jac2(i)*jac9(i))
271 jaci5=dett*(jac1(i)*jac9(i)-jac3(i)*jac7(i))
272 jaci8=dett*(jac2(i)*jac7(i)-jac1(i)*jac8(i))
273 jaci3=dett*(jac2(i)*jac6(i)-jac3(i)*jac5(i))
274 jaci6=dett*(jac3(i)*jac4(i)-jac1(i)*jac6(i))
275 jaci9=dett*(jac1(i)*jac5(i)-jac2(i)*jac4(i))
276
277 jac12=jaci1+jaci2
278 jac45=jaci4+jaci5
279 jac78=jaci7+jaci8
280
281 px1(i)=-jac12-jaci3
282 py1(i)=-jac45-jaci6
283 pz1(i)=-jac78-jaci9
284
285 px2(i)=-jac12+jaci3
286 py2(i)=-jac45+jaci6
287 pz2(i)=-jac78+jaci9
288
289 jac12=jaci1-jaci2
290 jac45=jaci4-jaci5
291 jac78=jaci7-jaci8
292
293 px3(i)= jac12+jaci3
294 py3(i)= jac45+jaci6
295 pz3(i)= jac78+jaci9
296 px4(i)= jac12-jaci3
297 py4(i)= jac45-jaci6
298 pz4(i)= jac78-jaci9
299
300 ENDDO
301
302 IF(jhbe/=0)THEN
303 DO i=1,nel
304
305 hx=(x1(i)-x2(i)+x3(i)-x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
306 hy=(y1(i)-y2(i)+y3(i)-y4(i)+y5(i)-y6(i)+y7
307 hz=(z1(i)-z2(i)+z3(i)-z4(i)+z5(i)-z6(i)+z7(i)-z8(i))
308 px1h1(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
309 px2h1(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
310 px3h1(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
311 px4h1(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
312 ENDDO
313
314 DO i=1,nel
315
316 hx=(x1(i)+x2(i)-x3(i)-x4(i)-x5(i)-x6(i)+x7(i)+x8(i))
317 hy=(y1(i)+y2(i)-y3(i)-y4(i)-y5(i)-y6(i)+y7(i)+y8(i))
318 hz=(z1(i)+z2(i)-z3(i)-z4(i)-z5(i)-z6(i)+z7(i)+z8(i))
319 px1h2(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
320 px2h2(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
321 px3h2(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
322 px4h2(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
323 ENDDO
324
325 DO i=1,nel
326
327 hx=(x1(i)-x2(i)-x3(i)+x4(i)-x5(i)+x6(i)+x7(i)-x8(i))
328 hy=(y1(i)-y2(i)-y3(i)+y4(i)-y5(i)+y6(i)+y7(i)-y8(i))
329 hz=(z1(i)-z2(i)-z3(i)+z4(i)-z5(i)+z6(i)+z7(i)-z8(i))
330 px1h3(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
331 px2h3(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
332 px3h3(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
333 px4h3(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
334 ENDDO
335 ENDIF
336
337 RETURN
subroutine schkjabt3(off, det, ngl, offg, nnega, index, nel, ismstr, jlag)